Hands-on lab — Day 100: Matrices and What They Represent

"""Exercise 1 — the from-scratch matrix, and the same thing in NumPy.

Run from the lab directory:

    .venv/bin/python3 examples/01_matrix_from_scratch.py

Every claim printed here is asserted immediately after it is printed, so the
script exits non-zero the moment a number stops being true.
"""

import numpy as np

from matrix import Matrix, ShapeMismatch

TOL = 1e-12  # every float comparison below uses numpy.allclose with this


def rule(title):
    print()
    print(title)
    print("-" * len(title))


A = Matrix([[2, 4, 1, 3], [0, 5, 2, 7], [6, 1, 4, 2]])
B = Matrix([[1, 0, 0, 1], [2, 2, 2, 2], [0, 3, 0, 3]])

npA = np.array(A.to_lists())
npB = np.array(B.to_lists())

rule("1a. Shape is (rows, columns) — rows first")
print(f"  from scratch : {A.shape}")
print(f"  numpy        : {npA.shape}")
print(f"  numpy ndim   : {npA.ndim}   (a matrix is a 2-dimensional array)")
print(f"  numpy size   : {npA.size}   (rows x columns, the entry count)")
assert A.shape == npA.shape == (3, 4)
assert npA.ndim == 2
assert npA.size == 12

rule("1b. Indexing counts from zero in both")
print("  the matrix:")
print(A.format())
print(f"  A[0, 0] = {A[0, 0]}   (row 0, column 0 — the top-left entry)")
print(f"  A[2, 3] = {A[2, 3]}   (row 2, column 3 — the bottom-right entry)")
print(f"  npA[2, 3] = {npA[2, 3]}")
print("  written in a paper, the bottom-right entry of this matrix is a_34,")
print("  because mathematics counts rows and columns from 1. Same entry.")
assert A[0, 0] == npA[0, 0] == 2
assert A[2, 3] == npA[2, 3] == 2
assert A.row(1) == list(npA[1]) == [0, 5, 2, 7]
assert A.col(1) == list(npA[:, 1]) == [4, 5, 1]

out_of_range = None
try:
    A[3, 0]
except IndexError as exc:
    out_of_range = str(exc)
print(f"  A[3, 0] raises IndexError: {out_of_range}")
assert out_of_range is not None and "shape (3, 4)" in out_of_range

rule("1c. Transpose swaps the axes")
print("  A.T:")
print(A.T.format())
print(f"  shape {A.shape} becomes {A.T.shape}")
assert A.T.shape == (4, 3)
assert np.array_equal(np.array(A.T.to_lists()), npA.T)
assert A.T.T == A, "transposing twice returns the original"

rule("1d. Addition is elementwise, and demands identical shapes")
print((A + B).format())
assert np.array_equal(np.array((A + B).to_lists()), npA + npB)

mismatch = None
try:
    A.add(Matrix([[1, 2], [3, 4]]))
except ShapeMismatch as exc:
    mismatch = str(exc)
print(f"  adding a (2, 2) raises ShapeMismatch: {mismatch}")
assert mismatch is not None
assert issubclass(ShapeMismatch, ValueError), "so `except ValueError` catches both"

rule("1e. Scalar multiplication multiplies every entry")
print((A * 3).format())
assert np.array_equal(np.array((A * 3).to_lists()), npA * 3)
assert np.array_equal(np.array((3 * A).to_lists()), 3 * npA)

rule("1f. The matrices worth recognising on sight")
print("  zeros(2, 3) — the additive nothing:")
print(Matrix.zeros(2, 3).format())
print("  identity(3) — leaves every vector exactly as it found it:")
print(Matrix.identity(3).format())
print("  diagonal([2, 5, 1]) — scales each coordinate by its own factor:")
print(Matrix.diagonal([2, 5, 1]).format())
S = Matrix([[1, 7, 3], [7, 4, 0], [3, 0, 9]])
print("  a symmetric matrix — equal to its own transpose:")
print(S.format())
print(f"  S.is_symmetric() = {S.is_symmetric()}   A.is_symmetric() = {A.is_symmetric()}")
assert np.array_equal(np.array(Matrix.identity(3).to_lists()), np.eye(3, dtype=int))
assert np.array_equal(np.array(Matrix.diagonal([2, 5, 1]).to_lists()), np.diag([2, 5, 1]))
assert S.is_symmetric() is True
assert A.is_symmetric() is False, "A is not even square, so it cannot be symmetric"

v = [1.0, 2.0, 3.0]
print(f"  identity(3) applied to {v} gives {Matrix.identity(3).apply_to(v)}")
assert np.allclose(Matrix.identity(3).apply_to(v), v, atol=TOL)

rule("1g. The one thing the from-scratch class cannot do")
print("  NumPy adds a (4,) row to every row of a (3, 4) array without asking:")
print(f"  npA + np.array([100, 200, 300, 400]) =\n{npA + np.array([100, 200, 300, 400])}")
refused = None
try:
    A.add([100, 200, 300, 400])
except TypeError as exc:
    refused = str(exc)
print(f"  the from-scratch class refuses: TypeError: {refused}")
assert refused is not None and "broadcasting" in refused
print("  That refusal is honest, and it is also the gap. Exercise 3 is about")
print("  what filling it costs you.")

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

"""Exercise 2 — one matrix, read three ways.

Run from the lab directory:

    .venv/bin/python3 examples/02_three_meanings.py

The same twelve numbers are read as a table of data, then as a collection of
vectors, then as a transformation. Nothing about the numbers changes. What
changes is which question you are asking, and every operation in this lab is
downstream of that choice.

The data is invented; see dataset.py.
"""

import numpy as np

from dataset import (
    COST_PER_BAG_PENCE,
    INGREDIENT_NAMES,
    LITRES_PER_BAG,
    LITRES_PER_INGREDIENT,
    MIX_NAMES,
    PRICE_PER_LITRE,
    RECIPES,
)
from matrix import Matrix

TOL = 1e-12


def L(a):
    """Render a numpy array as a plain Python list, with no dtype noise."""
    return np.asarray(a).tolist()


def rule(title):
    print()
    print(title)
    print("-" * len(title))


M = np.array(RECIPES)
prices = np.array(PRICE_PER_LITRE)

rule("Meaning 1 — a table of data: rows are items, columns are features")
header = "           " + "".join(f"{name:>9}" for name in INGREDIENT_NAMES)
print(header)
for name, row in zip(MIX_NAMES, M):
    print(f"{name:>10} " + "".join(f"{value:>9}" for value in row))
print()
print(f"  shape {M.shape}: {M.shape[0]} items described by {M.shape[1]} features.")
print("  This is the shape a CSV file arrives in (Day 65) and the shape a")
print("  SELECT returns (Day 85). Row 1 is a mix. Column 1 is a measurement")
print("  made of every mix. They are different kinds of thing living in one")
print("  rectangle, which is exactly why the axis argument keeps catching")
print("  people out.")
assert M.shape == (3, 4)
assert list(M[MIX_NAMES.index("Alpine")]) == [6, 1, 4, 2]
assert list(M[:, INGREDIENT_NAMES.index("grit")]) == [1, 2, 4]

rule("Meaning 2 — a collection of vectors, and WHICH collection matters")
print("  Read as rows, each vector is one mix in ingredient-space:")
for name, row in zip(MIX_NAMES, M):
    print(f"    {name:>10} -> {L(row)}   length {np.linalg.norm(row):.4f}")
print()
print("  Read as columns, each vector is one ingredient across the range:")
for name, col in zip(INGREDIENT_NAMES, M.T):
    print(f"    {name:>10} -> {L(col)}   length {np.linalg.norm(col):.4f}")
print()
print("  Same twelve numbers. Three vectors of length 4, or four vectors of")
print("  length 3 — and the norms are entirely different numbers, so an")
print("  operation that assumed the wrong one gives a plausible wrong answer")
print("  rather than an error.")
row_norms = np.linalg.norm(M, axis=1)
col_norms = np.linalg.norm(M, axis=0)
assert row_norms.shape == (3,)
assert col_norms.shape == (4,)
# Seedling: sqrt(4 + 16 + 1 + 9) = sqrt(30); worked out by hand.
assert np.allclose(row_norms[0], np.sqrt(30.0), atol=TOL)
# base column: sqrt(4 + 0 + 36) = sqrt(40)
assert np.allclose(col_norms[0], np.sqrt(40.0), atol=TOL)
assert not np.allclose(np.sort(row_norms), np.sort(col_norms[:3]), atol=TOL)

rule("Meaning 3 — a transformation: give it a vector, get a vector back")
print(f"  in : the price of each ingredient in pence per litre {L(prices)}")
print("       (a vector of length 4, one entry per COLUMN of the matrix)")
by_hand = []
for name, row in zip(MIX_NAMES, M):
    terms = " + ".join(f"{q}*{p}" for q, p in zip(row, prices))
    total = int(np.sum(row * prices))
    by_hand.append(total)
    print(f"    {name:>10}: {terms} = {total}")
print(f"  out: the cost of one bag of each mix in pence {by_hand}")
print("       (a vector of length 3, one entry per ROW of the matrix)")
print()
print("  A (3, 4) matrix eats a vector of length 4 and returns one of length 3.")
print("  The 4 is consumed; the 3 survives. That is the whole shape rule, and")
print("  it is why a shape error is always a disagreement about which of the")
print("  three meanings each side assumed.")
assert by_hand == COST_PER_BAG_PENCE

scratch = Matrix(RECIPES).apply_to(PRICE_PER_LITRE)
broadcast = (M * prices).sum(axis=1)
print()
print(f"  from-scratch double loop : {scratch}")
print(f"  numpy (M * prices).sum(axis=1): {L(broadcast)}")
print("  Both are the same arithmetic. The packed operator that names it,")
print("  M @ prices, is Day 101.")
assert scratch == COST_PER_BAG_PENCE
assert np.allclose(broadcast, COST_PER_BAG_PENCE, atol=TOL)

rule("The identity transformation, on this data")
identity = np.eye(4, dtype=int)
print("  A (4, 4) identity matrix applied to the price vector returns it:")
print(f"    {L((identity * prices).sum(axis=1))}")
assert np.allclose((identity * prices).sum(axis=1), prices, atol=TOL)
print("  A diagonal matrix instead scales each ingredient's price separately.")
doubler = np.diag([2, 1, 1, 1])
print(f"    doubling only the base price: {L((doubler * prices).sum(axis=1))}")
assert list((doubler * prices).sum(axis=1)) == [20, 2, 5, 1]

rule("The two totals the table meaning wants")
print(f"  litres in each bag (sum across each row)      : {L(M.sum(axis=1))}")
print(f"  litres of each ingredient (sum down each col) : {L(M.sum(axis=0))}")
assert list(M.sum(axis=1)) == LITRES_PER_BAG
assert list(M.sum(axis=0)) == LITRES_PER_INGREDIENT

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

"""Exercise 3 — reshape, flatten, slice, and the word "view".

Run from the lab directory:

    .venv/bin/python3 examples/03_views_and_copies.py

A NumPy array is two things: a flat block of memory, and a description of how
to read it as a grid. Reshaping usually rewrites only the description. When it
does, the two names are two windows onto ONE block of numbers, and writing
through either window changes what the other one sees.

That is not a bug and it is not a subtlety you can ignore. It is the reason
NumPy is fast, and it is the reason a function that "just reshapes" its input
can quietly edit your dataset.
"""

import numpy as np

from dataset import RECIPES

TOL = 1e-12


def L(a):
    """Render a numpy array as a plain Python list, with no dtype noise."""
    return np.asarray(a).tolist()


def rule(title):
    print()
    print(title)
    print("-" * len(title))


def fresh():
    return np.array(RECIPES)


rule("3a. Reshape rewrites the description, not the numbers")
M = fresh()
flat = M.reshape(12)
print(f"  M.shape    = {M.shape}")
print(f"  flat.shape = {flat.shape}")
print(f"  flat       = {L(flat)}")
print("  The entries come out row by row: the whole of row 0, then row 1, then")
print("  row 2. NumPy calls that C order, and it is the default everywhere.")
assert list(flat) == [2, 4, 1, 3, 0, 5, 2, 7, 6, 1, 4, 2]
assert M.size == flat.size == 12

print()
print(f"  M.reshape(2, 6) =\n{M.reshape(2, 6)}")
print(f"  M.reshape(4, 3) =\n{M.reshape(4, 3)}")
print("  Any shape whose entries multiply to 12 is allowed, and -1 means")
print(f"  'work it out': M.reshape(6, -1).shape = {M.reshape(6, -1).shape}")
assert M.reshape(6, -1).shape == (6, 2)

bad = None
try:
    M.reshape(5, 3)
except ValueError as exc:
    bad = str(exc)
print(f"  M.reshape(5, 3) raises ValueError: {bad}")
assert bad is not None and "reshape" in bad

rule("3b. The proof: mutate the reshaped array, watch the original change")
M = fresh()
flat = M.reshape(12)
print(f"  before: M[0, 0] = {M[0, 0]}, flat[0] = {flat[0]}")
flat[0] = 99
print("  flat[0] = 99")
print(f"  after : M[0, 0] = {M[0, 0]}, flat[0] = {flat[0]}")
print("  Nothing was assigned to M. M changed anyway.")
assert M[0, 0] == 99
print(f"  flat.base is M           -> {flat.base is M}")
print(f"  numpy.shares_memory(M, flat) -> {np.shares_memory(M, flat)}")
assert flat.base is M
assert np.shares_memory(M, flat)

rule("3c. .copy() breaks the link")
M = fresh()
independent = M.copy().reshape(12)
print(f"  before: M[0, 0] = {M[0, 0]}, independent[0] = {independent[0]}")
independent[0] = 99
print("  independent[0] = 99")
print(f"  after : M[0, 0] = {M[0, 0]}, independent[0] = {independent[0]}")
assert M[0, 0] == 2, "the original is untouched"
assert independent[0] == 99
print(f"  numpy.shares_memory(M, independent) -> {np.shares_memory(M, independent)}")
assert not np.shares_memory(M, independent)
print()
print("  A warning about the obvious-looking test. `.base is None` is NOT the")
print("  question to ask:")
print(f"    independent.base is None -> {independent.base is None}")
print("  It is False, and the array is still fully independent of M. The reason")
print("  is that M.copy().reshape(12) made TWO arrays: an anonymous copy, and a")
print("  reshaped view of that copy. `.base` truthfully points at the copy,")
print("  which no longer has a name and which nothing else can reach.")
print("  numpy.shares_memory(a, b) asks the question you actually care about.")
assert independent.base is not None
assert np.shares_memory(independent, independent.base)
assert not np.shares_memory(M, independent.base)
print("  A view costs nothing and shares everything. A copy costs the memory")
print("  and shares nothing. Neither is the right answer; knowing which one")
print("  you have is.")

rule("3d. A slice is a view too — this is the one that bites")
M = fresh()
grit = M[:, 2]
print(f"  grit = M[:, 2] = {L(grit)}   (the third column, counting from 0)")
grit[0] = 50
print("  grit[0] = 50")
print(f"  M is now:\n{M}")
assert M[0, 2] == 50
assert np.shares_memory(M, grit)
print("  Slicing a Python list gives you a new list. Slicing a NumPy array")
print("  gives you a window. The syntax is identical and the behaviour is not.")

py_list = [[2, 4, 1, 3], [0, 5, 2, 7]]
sliced = py_list[0][:]
sliced[0] = 50
print(f"  for contrast, a Python list slice: original still {py_list[0]}")
assert py_list[0][0] == 2

rule("3e. Fancy indexing always copies")
M = fresh()
picked = M[[0, 2]]
picked[0, 0] = 77
print("  M[[0, 2]] selects rows 0 and 2 by a LIST of indices, not a slice.")
print(f"  after picked[0, 0] = 77, M[0, 0] is still {M[0, 0]}")
assert M[0, 0] == 2
assert not np.shares_memory(M, picked)
print(f"  numpy.shares_memory(M, picked) -> {np.shares_memory(M, picked)}")
print("  Rule of thumb: basic slicing with colons gives a view; indexing with")
print("  a list or a boolean mask gives a copy. When in doubt, ask")
print("  numpy.shares_memory rather than guessing.")

rule("3f. .ravel() versus .flatten()")
M = fresh()
r = M.ravel()
f = M.flatten()
print(f"  M.ravel()   shares memory with M: {np.shares_memory(M, r)}")
print(f"  M.flatten() shares memory with M: {np.shares_memory(M, f)}")
print("  Same numbers, same shape, opposite ownership. ravel gives a view when")
print("  it can; flatten always copies. The names do not tell you that, so this")
print("  is a fact to memorise rather than derive.")
assert np.shares_memory(M, r)
assert not np.shares_memory(M, f)
assert np.array_equal(r, f)

rule("3g. Transpose is a view as well")
M = fresh()
t = M.T
print(f"  M.T.shape = {t.shape}, shares memory with M: {np.shares_memory(M, t)}")
t[0, 0] = 42
print(f"  after t[0, 0] = 42, M[0, 0] = {M[0, 0]}")
assert np.shares_memory(M, t)
assert M[0, 0] == 42
print("  Transposing a large matrix moves no numbers at all — it swaps two")
print("  entries in the description of how to walk the memory. That is why")
print("  transpose is free and why it is a view.")

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

"""Exercise 4 — broadcasting: the rule, a success, a failure, and a trap.

Run from the lab directory:

    .venv/bin/python3 examples/04_broadcasting.py

Broadcasting is where NumPy stops being obvious. Up to here, every operation
did the boring thing. Broadcasting invents entries that were never written
down, and it does so silently, which makes it the single most productive
source of results that are wrong without being errors.

The rule, from the NumPy documentation, applied right to left:

    1. Line the two shapes up from the RIGHT-hand end.
    2. A missing entry on the left of the shorter shape counts as 1.
    3. Two dimensions are compatible when they are equal, or one of them is 1.
    4. If any pair is neither, the operation is a ValueError.
    5. The result takes the larger of each pair.

Nothing is copied. The stretching is a fiction maintained by walking the
smaller array's memory more than once.
"""

import numpy as np

from dataset import PRICE_PER_LITRE, RECIPES

TOL = 1e-12


def L(a):
    """Render a numpy array as a plain Python list, with no dtype noise."""
    return np.asarray(a).tolist()


def rule(title):
    print()
    print(title)
    print("-" * len(title))


def align(a, b):
    """Print the right-aligned shape comparison the rule actually describes."""
    width = max(len(a), len(b))
    pad = lambda s: (1,) * (width - len(s)) + s  # noqa: E731 - a local, read once
    top, bottom = pad(a), pad(b)
    ok = all(x == y or x == 1 or y == 1 for x, y in zip(top, bottom))
    result = tuple(max(x, y) for x, y in zip(top, bottom)) if ok else None
    print(f"    {str(a):>12}   padded to {top}")
    print(f"    {str(b):>12}   padded to {bottom}")
    print(f"    -> {'result ' + str(result) if ok else 'INCOMPATIBLE'}")
    return result


M = np.array(RECIPES)
prices = np.array(PRICE_PER_LITRE)

rule("4a. The simplest case: an array and one number")
print(f"  M * 2 has shape {(M * 2).shape}; the 2 is treated as (1, 1)")
print(M * 2)
assert np.array_equal(M * 2, M + M)

rule("4b. (3, 4) with (4,) — this succeeds, and here is why")
align((3, 4), (4,))
print()
print(f"  M       shape {M.shape}:\n{M}")
print(f"  prices  shape {prices.shape}: {L(prices)}")
scaled = M * prices
print(f"  M * prices shape {scaled.shape}:\n{scaled}")
print()
print("  Read row by row: the price vector was applied to EVERY row, because")
print("  the rows are the axis that had to stretch. Row 0 became")
print(f"  {L(M[0])} * {L(prices)} = {L(scaled[0])}.")
assert scaled.shape == (3, 4)
assert list(scaled[0]) == [20, 8, 5, 3]
assert list(scaled[1]) == [0, 10, 10, 7]
assert list(scaled[2]) == [60, 2, 20, 2]
print()
print("  And nothing was copied. numpy.broadcast_to shows the fiction directly:")
stretched = np.broadcast_to(prices, (3, 4))
print(stretched)
print(f"  shares memory with prices: {np.shares_memory(stretched, prices)}")
print(f"  writeable: {stretched.flags.writeable}  (it has to be read-only —")
print("  one write would appear in three places at once)")
assert np.shares_memory(stretched, prices)
assert stretched.flags.writeable is False

rule("4c. (3, 4) with (3,) — this fails, and the rule predicted it")
align((3, 4), (3,))
print()
per_mix = np.array([100, 200, 300])
failure = None
try:
    M + per_mix
except ValueError as exc:
    failure = exc
print(f"  M + numpy.array({L(per_mix)}) raises")
print(f"    {type(failure).__name__}: {failure}")
assert isinstance(failure, ValueError)
assert "could not be broadcast together" in str(failure)
print()
print("  There is nothing special about 3 versus 4 here. The trailing")
print("  dimensions were 4 and 3, neither equal nor 1, so the rule stopped.")
print("  The fix is to say which axis you meant:")
column = per_mix.reshape(3, 1)
align((3, 4), (3, 1))
print(f"  per_mix.reshape(3, 1) has shape {column.shape}, and M + it works:")
print(M + column)
assert (M + column).shape == (3, 4)
assert list((M + column)[0]) == [102, 104, 101, 103]
assert list((M + column)[2]) == [306, 301, 304, 302]

rule("4d. The trap: a square matrix, where the wrong answer is not an error")
S = np.array(
    [
        [1.0, 2.0, 3.0, 4.0],
        [10.0, 20.0, 30.0, 40.0],
        [100.0, 200.0, 300.0, 400.0],
        [1000.0, 2000.0, 3000.0, 4000.0],
    ]
)
print("  A (4, 4) matrix. Say you want to centre each ROW on its own mean —")
print("  a completely ordinary preprocessing step.")
row_means = S.mean(axis=1)
print(f"  row means, shape {row_means.shape}: {L(row_means)}")
wrong = S - row_means
right = S - S.mean(axis=1, keepdims=True)
print()
print(f"  S - row_means            (shape {wrong.shape}) — no error, no warning:")
print(wrong)
print(f"  S - S.mean(axis=1, keepdims=True) (shape {right.shape}) — what you meant:")
print(right)
print()
print("  The first one subtracted row 0's mean from COLUMN 0, row 1's mean from")
print("  column 1, and so on, because a (4,) lines up against the last axis and")
print("  the last axis is columns. It is a transposed answer wearing the right")
print("  shape. On a non-square matrix it would have been a ValueError and you")
print("  would have found out in one second.")
assert not np.allclose(wrong, right, atol=TOL)
# The correct centring makes every row sum to zero. The wrong one does not.
assert np.allclose(right.sum(axis=1), np.zeros(4), atol=1e-9)
assert not np.allclose(wrong.sum(axis=1), np.zeros(4), atol=1e-9)
print("  The check that catches it: after centring rows, every ROW must sum to 0.")
print(f"    right.sum(axis=1) = {L(right.sum(axis=1))}")
print(f"    wrong.sum(axis=1) = {L(wrong.sum(axis=1))}")

rule("4e. keepdims is the habit that prevents the whole family of bugs")
print(f"  S.mean(axis=1).shape                 = {S.mean(axis=1).shape}")
print(f"  S.mean(axis=1, keepdims=True).shape  = {S.mean(axis=1, keepdims=True).shape}")
print(f"  S.mean(axis=0, keepdims=True).shape  = {S.mean(axis=0, keepdims=True).shape}")
print("  keepdims=True leaves a 1 where the axis was, so the result still says")
print("  out loud which axis it came from, and broadcasting lines it up the way")
print("  you intended instead of the way the padding rule happened to choose.")
assert S.mean(axis=1, keepdims=True).shape == (4, 1)
assert S.mean(axis=0, keepdims=True).shape == (1, 4)

rule("4f. The outer-product surprise: (3, 1) with (1, 4)")
a = np.array([1, 2, 3]).reshape(3, 1)
b = np.array([10, 20, 30, 40]).reshape(1, 4)
align((3, 1), (1, 4))
print(f"  a * b has shape {(a * b).shape} — twelve entries from seven numbers:")
print(a * b)
print("  Both sides stretched. If you expected a length-3 or length-4 answer,")
print("  the size of the output is the first thing that tells you otherwise.")
assert (a * b).shape == (3, 4)
assert list((a * b)[2]) == [30, 60, 90, 120]

rule("4g. How to check your shapes before you trust the numbers")
print("  numpy.broadcast_shapes answers the question without doing the work:")
print(f"    broadcast_shapes((3, 4), (4,))    = {np.broadcast_shapes((3, 4), (4,))}")
print(f"    broadcast_shapes((3, 4), (3, 1))  = {np.broadcast_shapes((3, 4), (3, 1))}")
predicted = None
try:
    np.broadcast_shapes((3, 4), (3,))
except ValueError as exc:
    predicted = str(exc)
print(f"    broadcast_shapes((3, 4), (3,))    raises ValueError: {predicted}")
assert np.broadcast_shapes((3, 4), (4,)) == (3, 4)
assert predicted is not None
print("  Print .shape at every step you are unsure about. It costs one line and")
print("  it is the only debugging tool that works before the numbers are wrong.")

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

"""Exercise 5 — settling axis=0 against axis=1, once, with a small example.

Run from the lab directory:

    .venv/bin/python3 examples/05_axes.py

Everybody gets this wrong at least once, and the reason is that both readings
sound right in English. "Sum along the rows" can mean sum each row, or sum in
the direction the rows are stacked. English will not save you. One rule will:

    THE AXIS YOU NAME IS THE AXIS THAT DISAPPEARS.

A (3, 4) array summed with axis=0 loses the 3 and returns shape (4,).
A (3, 4) array summed with axis=1 loses the 4 and returns shape (3,).

Because 3 and 4 are different here, the shape of the answer tells you which
one you got, every time, without thinking. That is why this matrix is not
square.
"""

import numpy as np

from dataset import (
    INGREDIENT_NAMES,
    LITRES_PER_BAG,
    LITRES_PER_INGREDIENT,
    MIX_NAMES,
    RECIPES,
)

TOL = 1e-12


def L(a):
    """Render a numpy array as a plain Python list, with no dtype noise."""
    return np.asarray(a).tolist()


def rule(title):
    print()
    print(title)
    print("-" * len(title))


M = np.array(RECIPES)

rule("5a. The matrix, small enough to check on paper")
print("           " + "".join(f"{n:>9}" for n in INGREDIENT_NAMES))
for name, row in zip(MIX_NAMES, M):
    print(f"{name:>10} " + "".join(f"{v:>9}" for v in row))
print(f"  shape {M.shape}  -> axis 0 has length 3, axis 1 has length 4")

rule("5b. axis=0 removes axis 0, which is the rows")
total0 = M.sum(axis=0)
print(f"  M.sum(axis=0) = {L(total0)}   shape {total0.shape}")
print("  Four numbers, one per COLUMN. The three mixes were collapsed into one.")
print("  Worked by hand, column by column:")
for name, col in zip(INGREDIENT_NAMES, M.T):
    print(f"    {name:>10}: {' + '.join(str(v) for v in col)} = {sum(col)}")
assert total0.shape == (4,)
assert list(total0) == LITRES_PER_INGREDIENT == [8, 10, 7, 12]

rule("5c. axis=1 removes axis 1, which is the columns")
total1 = M.sum(axis=1)
print(f"  M.sum(axis=1) = {L(total1)}   shape {total1.shape}")
print("  Three numbers, one per ROW. The four ingredients were collapsed into one.")
print("  Worked by hand, row by row:")
for name, row in zip(MIX_NAMES, M):
    print(f"    {name:>10}: {' + '.join(str(v) for v in row)} = {sum(row)}")
assert total1.shape == (3,)
assert list(total1) == LITRES_PER_BAG == [10, 14, 13]

rule("5d. No axis at all collapses everything")
print(f"  M.sum() = {M.sum()}   shape {np.shape(M.sum())}  (a scalar)")
print("  Which is also the sum of either of the two answers above:")
print(f"    sum of the axis=0 answer: {total0.sum()}")
print(f"    sum of the axis=1 answer: {total1.sum()}")
assert M.sum() == 37
assert total0.sum() == total1.sum() == M.sum()

rule("5e. The same rule holds for every reduction, not just sum")
table = [
    ("sum", np.sum),
    ("mean", np.mean),
    ("min", np.min),
    ("max", np.max),
    ("argmax", np.argmax),
]
print(f"  {'function':>10} {'axis=0 result':>34} {'shape':>8}")
for name, fn in table:
    value = fn(M, axis=0)
    print(f"  {name:>10} {str(np.round(value, 4)):>34} {str(value.shape):>8}")
print(f"  {'function':>10} {'axis=1 result':>34} {'shape':>8}")
for name, fn in table:
    value = fn(M, axis=1)
    print(f"  {name:>10} {str(np.round(value, 4)):>34} {str(value.shape):>8}")
assert np.allclose(M.mean(axis=0), [8 / 3, 10 / 3, 7 / 3, 4.0], atol=TOL)
assert np.allclose(M.mean(axis=1), [2.5, 3.5, 3.25], atol=TOL)
assert list(np.argmax(M, axis=0)) == [2, 1, 2, 1]
assert list(np.argmax(M, axis=1)) == [1, 3, 0]

rule("5f. argmax with an axis returns positions, not values")
print(f"  numpy.argmax(M, axis=1) = {L(np.argmax(M, axis=1))}")
print("  Read it as: within row 0 the largest entry is at column 1; within")
print("  row 1 at column 3; within row 2 at column 0. Turning those into names:")
for name, j in zip(MIX_NAMES, np.argmax(M, axis=1)):
    print(f"    {name:>10}'s largest ingredient is {INGREDIENT_NAMES[j]} ({M[MIX_NAMES.index(name), j]} litres)")
assert INGREDIENT_NAMES[np.argmax(M, axis=1)[2]] == "base"

rule("5g. The same rule tells you what keepdims does")
print(f"  M.sum(axis=0).shape                = {M.sum(axis=0).shape}")
print(f"  M.sum(axis=0, keepdims=True).shape = {M.sum(axis=0, keepdims=True).shape}")
print(f"  M.sum(axis=1, keepdims=True).shape = {M.sum(axis=1, keepdims=True).shape}")
print("  keepdims=True leaves a 1 in place of the axis instead of removing it,")
print("  which is exactly what broadcasting needs to line the answer back up")
print("  against the matrix it came from.")
assert M.sum(axis=0, keepdims=True).shape == (1, 4)
assert M.sum(axis=1, keepdims=True).shape == (3, 1)
share = M / M.sum(axis=1, keepdims=True)
print("  Each row as a fraction of its own total — the operation that")
print("  keepdims exists for:")
print(np.round(share, 4))
assert np.allclose(share.sum(axis=1), np.ones(3), atol=TOL)

rule("5h. The table worth memorising")
print("  For a 2-D array of shape (rows, columns):")
print()
print("    argument   collapses   answer is one number per   answer shape")
print("    axis=0     the rows    column                     (columns,)")
print("    axis=1     the columns row                        (rows,)")
print("    no axis    everything  whole array                ()")
print()
print("  And the check that never lies: look at the LENGTH of the answer.")
print(f"  Here, a length-{len(total0)} answer came from axis=0 and a length-{len(total1)}")
print("  answer came from axis=1, because 4 columns and 3 rows.")

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

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

Both `examples/` and `starter/` contain a module called `matrix`, 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
`matrix` 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
`matrix` or `dataset` 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 ("matrix", "dataset"):
    module = sys.modules.get(name)
    origin = getattr(module, "__file__", "") or ""
    if module is not None and not origin.startswith(HERE):
        del sys.modules[name]

"""The one small dataset this lab reads three different ways.

INVENTED DATA. Fenwick Road Garden Centre does not exist, and neither do
these recipes or these prices. The numbers were chosen so that every result
in this lab can be checked by hand on paper in under a minute, which is the
only property that matters for a teaching example.

Three potting mixes, four ingredients, litres of each ingredient per bag:

                base   bark   grit   compost
    Seedling       2      4      1         3
    Container      0      5      2         7
    Alpine         6      1      4         2

and the ingredient prices, in pence per litre:

    base 10, bark 2, grit 5, compost 1

Everything else in the lab is derived from those two blocks of numbers.
"""

MIX_NAMES = ["Seedling", "Container", "Alpine"]
INGREDIENT_NAMES = ["base", "bark", "grit", "compost"]

# The recipes, as a plain nested list: one inner list per mix.
RECIPES = [
    [2, 4, 1, 3],
    [0, 5, 2, 7],
    [6, 1, 4, 2],
]

# Pence per litre, one entry per ingredient, in the same column order.
PRICE_PER_LITRE = [10, 2, 5, 1]

# The answers, worked out by hand, so the code has something to be checked
# against rather than merely something to print.
#
#   Seedling : 2*10 + 4*2 + 1*5 + 3*1 = 20 +  8 +  5 + 3 = 36
#   Container: 0*10 + 5*2 + 2*5 + 7*1 =  0 + 10 + 10 + 7 = 27
#   Alpine   : 6*10 + 1*2 + 4*5 + 2*1 = 60 +  2 + 20 + 2 = 84
COST_PER_BAG_PENCE = [36, 27, 84]

# Row sums: how many litres are in each bag.
LITRES_PER_BAG = [10, 14, 13]

# Column sums: how many litres of each ingredient one bag of every mix needs.
LITRES_PER_INGREDIENT = [8, 10, 7, 12]

"""A matrix built from first principles, on nothing but nested lists.

This is the reference implementation the lab's tests assert against, and the
thing NumPy is compared to. It deliberately implements only five ideas:

    * shape, as a (rows, columns) tuple;
    * indexing with a (row, column) pair, counting from zero;
    * transpose, which swaps the two axes;
    * addition, which is elementwise and requires identical shapes;
    * scalar multiplication, which multiplies every entry by one number.

It deliberately does NOT implement broadcasting. That omission is the point
of the last exercise: broadcasting is the first thing on this list that is
genuinely difficult to write and genuinely easy to misuse, and seeing the
gap is how you learn what NumPy is doing for you.

Every operation returns a NEW Matrix. Nothing here mutates in place, which
makes the class boring and safe, and makes the contrast with NumPy views
(where mutation travels) sharper in exercise 4.
"""

from __future__ import annotations


class ShapeMismatch(ValueError):
    """Raised when two matrices cannot be combined because their shapes differ.

    A subclass of ValueError on purpose: NumPy raises ValueError for the same
    situation, so code that catches ValueError catches both. The tests assert
    that relationship rather than assuming it.
    """


class Matrix:
    """A rectangular grid of numbers, stored as a list of row lists."""

    def __init__(self, rows):
        rows = [list(row) for row in rows]
        if not rows:
            raise ValueError("a matrix needs at least one row")
        width = len(rows[0])
        if width == 0:
            raise ValueError("a matrix needs at least one column")
        for i, row in enumerate(rows):
            if len(row) != width:
                raise ValueError(
                    f"row {i} has {len(row)} entries but row 0 has {width}; "
                    "a matrix is rectangular, so every row must be the same length"
                )
        self._rows = rows

    # -- shape ------------------------------------------------------------

    @property
    def shape(self):
        """(rows, columns) — rows first, because that is the convention."""
        return (len(self._rows), len(self._rows[0]))

    @property
    def n_rows(self):
        return self.shape[0]

    @property
    def n_cols(self):
        return self.shape[1]

    # -- indexing ---------------------------------------------------------

    def __getitem__(self, position):
        """m[i, j] — row i, column j, both counting from 0."""
        i, j = self._check_position(position)
        return self._rows[i][j]

    def __setitem__(self, position, value):
        i, j = self._check_position(position)
        self._rows[i][j] = value

    def _check_position(self, position):
        if not isinstance(position, tuple) or len(position) != 2:
            raise TypeError(
                "index a Matrix with a (row, column) pair, for example m[0, 2]"
            )
        i, j = position
        n_rows, n_cols = self.shape
        if not (0 <= i < n_rows) or not (0 <= j < n_cols):
            raise IndexError(
                f"position {position} is outside a matrix of shape {self.shape}; "
                f"valid rows are 0..{n_rows - 1} and valid columns are 0..{n_cols - 1}"
            )
        return i, j

    def row(self, i):
        """Row i as a plain list — a copy, so editing it changes nothing here."""
        return list(self._rows[i])

    def col(self, j):
        """Column j as a plain list. Note this has to walk every row."""
        return [row[j] for row in self._rows]

    # -- operations -------------------------------------------------------

    def transpose(self):
        """Swap rows and columns: an (r, c) matrix becomes (c, r)."""
        n_rows, n_cols = self.shape
        return Matrix([[self._rows[i][j] for i in range(n_rows)] for j in range(n_cols)])

    @property
    def T(self):
        """Spelled the way NumPy spells it, so the comparison reads cleanly."""
        return self.transpose()

    def add(self, other):
        """Elementwise addition. Both shapes must match exactly — no broadcasting."""
        if not isinstance(other, Matrix):
            raise TypeError(
                "add expects another Matrix; this class has no broadcasting, "
                "so a plain number or a list is not accepted"
            )
        if self.shape != other.shape:
            raise ShapeMismatch(
                f"cannot add {self.shape} to {other.shape}: "
                "addition here is entry by entry, so the shapes must be identical"
            )
        n_rows, n_cols = self.shape
        return Matrix(
            [[self[i, j] + other[i, j] for j in range(n_cols)] for i in range(n_rows)]
        )

    def __add__(self, other):
        return self.add(other)

    def scale(self, k):
        """Multiply every entry by the number k."""
        n_rows, n_cols = self.shape
        return Matrix([[self[i, j] * k for j in range(n_cols)] for i in range(n_rows)])

    def __mul__(self, k):
        return self.scale(k)

    __rmul__ = __mul__

    def apply_to(self, vector):
        """Treat the matrix as a transformation and apply it to one vector.

        Each output entry is the sum of one row multiplied entry by entry
        against the input vector. An (r, c) matrix therefore eats a vector of
        length c and returns a vector of length r.

        This is written out as an explicit double loop on purpose: the packed
        name for it, and the operator NumPy spells `@`, is Day 101's subject.
        """
        vector = list(vector)
        n_rows, n_cols = self.shape
        if len(vector) != n_cols:
            raise ShapeMismatch(
                f"a {self.shape} matrix transforms a vector of length {n_cols}, "
                f"but this vector has length {len(vector)}"
            )
        return [
            sum(self[i, j] * vector[j] for j in range(n_cols)) for i in range(n_rows)
        ]

    # -- constructors for the matrices worth recognising on sight ---------

    @classmethod
    def zeros(cls, n_rows, n_cols):
        return cls([[0] * n_cols for _ in range(n_rows)])

    @classmethod
    def identity(cls, n):
        """The matrix that leaves every vector exactly as it found it."""
        return cls([[1 if i == j else 0 for j in range(n)] for i in range(n)])

    @classmethod
    def diagonal(cls, values):
        values = list(values)
        n = len(values)
        return cls([[values[i] if i == j else 0 for j in range(n)] for i in range(n)])

    # -- reporting --------------------------------------------------------

    def is_symmetric(self):
        """True when the matrix equals its own transpose. Requires squareness."""
        n_rows, n_cols = self.shape
        if n_rows != n_cols:
            return False
        return all(
            self[i, j] == self[j, i] for i in range(n_rows) for j in range(n_cols)
        )

    def to_lists(self):
        """A plain nested list — the shape numpy.array() wants."""
        return [list(row) for row in self._rows]

    def __eq__(self, other):
        if not isinstance(other, Matrix):
            return NotImplemented
        return self._rows == other._rows

    def __repr__(self):
        body = ", ".join(repr(row) for row in self._rows)
        return f"Matrix([{body}])  # shape {self.shape}"

    def format(self, width=4):
        """A readable grid, one row per line, right-aligned in fixed columns."""
        return "\n".join(
            " ".join(f"{value:>{width}}" for value in row) for row in self._rows
        )

"""The reference test suite: real values, real shapes, one stated tolerance.

Run from the lab directory:

    .venv/bin/pytest examples

Every float comparison here goes through numpy.allclose with atol=TOL below.
Integer and shape comparisons are exact, because there is nothing to round.
"""

import numpy as np
import pytest

from dataset import (
    COST_PER_BAG_PENCE,
    INGREDIENT_NAMES,
    LITRES_PER_BAG,
    LITRES_PER_INGREDIENT,
    MIX_NAMES,
    PRICE_PER_LITRE,
    RECIPES,
)
from matrix import Matrix, ShapeMismatch

# The stated tolerance. 1e-12 is far tighter than any error these numbers can
# accumulate — they are small integers and exact halves — and stating it is
# the point: a float comparison without a declared tolerance is a guess.
TOL = 1e-12


@pytest.fixture
def M():
    """A fresh (3, 4) array every test, because several tests mutate views."""
    return np.array(RECIPES)


@pytest.fixture
def A():
    return Matrix(RECIPES)


# --------------------------------------------------------------------------
# 1. The from-scratch class, asserted against NumPy
# --------------------------------------------------------------------------


def test_shape_is_rows_then_columns(A, M):
    assert A.shape == (3, 4)
    assert A.shape == M.shape
    assert A.n_rows == 3 and A.n_cols == 4
    assert M.ndim == 2 and M.size == 12


def test_indexing_counts_from_zero(A, M):
    assert A[0, 0] == 2 and M[0, 0] == 2
    assert A[2, 3] == 2 and M[2, 3] == 2
    assert A[1, 3] == 7 and M[1, 3] == 7


def test_indexing_out_of_range_is_an_IndexError(A):
    with pytest.raises(IndexError) as caught:
        A[3, 0]
    assert "shape (3, 4)" in str(caught.value)


def test_indexing_with_a_single_number_is_rejected(A):
    with pytest.raises(TypeError):
        A[0]


def test_rows_and_columns_agree_with_numpy(A, M):
    assert A.row(1) == M[1].tolist() == [0, 5, 2, 7]
    assert A.col(1) == M[:, 1].tolist() == [4, 5, 1]


def test_transpose_swaps_the_axes(A, M):
    assert A.T.shape == (4, 3) == M.T.shape
    assert A.T.to_lists() == M.T.tolist()
    assert A.T.T == A


def test_addition_is_elementwise(A, M):
    B = Matrix([[1, 0, 0, 1], [2, 2, 2, 2], [0, 3, 0, 3]])
    npB = np.array(B.to_lists())
    assert (A + B).to_lists() == (M + npB).tolist()
    assert (A + B)[0, 0] == 3


def test_addition_of_different_shapes_raises(A):
    with pytest.raises(ShapeMismatch):
        A.add(Matrix([[1, 2], [3, 4]]))
    # And, because ShapeMismatch subclasses ValueError, this catches it too —
    # which is the same exception type NumPy raises for the same mistake.
    with pytest.raises(ValueError):
        A.add(Matrix([[1, 2], [3, 4]]))


def test_scalar_multiplication_from_both_sides(A, M):
    assert (A * 3).to_lists() == (M * 3).tolist()
    assert (3 * A).to_lists() == (3 * M).tolist()
    assert (A * 0).to_lists() == Matrix.zeros(3, 4).to_lists()


def test_identity_and_diagonal_match_numpy():
    assert Matrix.identity(3).to_lists() == np.eye(3, dtype=int).tolist()
    assert Matrix.diagonal([2, 5, 1]).to_lists() == np.diag([2, 5, 1]).tolist()


def test_symmetry(A):
    S = Matrix([[1, 7, 3], [7, 4, 0], [3, 0, 9]])
    assert S.is_symmetric()
    assert S.to_lists() == S.T.to_lists()
    assert not A.is_symmetric(), "a (3, 4) matrix is not square, so not symmetric"


def test_identity_leaves_a_vector_alone():
    v = [1.5, -2.0, 0.25]
    assert np.allclose(Matrix.identity(3).apply_to(v), v, atol=TOL)


def test_from_scratch_class_cannot_broadcast(A):
    with pytest.raises(TypeError):
        A.add([100, 200, 300, 400])


# --------------------------------------------------------------------------
# 2. One matrix, three meanings
# --------------------------------------------------------------------------


def test_meaning_one_table_rows_are_items_columns_are_features(M):
    assert M.shape == (len(MIX_NAMES), len(INGREDIENT_NAMES))
    assert M[MIX_NAMES.index("Alpine")].tolist() == [6, 1, 4, 2]
    assert M[:, INGREDIENT_NAMES.index("grit")].tolist() == [1, 2, 4]


def test_meaning_two_rows_and_columns_are_different_vector_sets(M):
    row_norms = np.linalg.norm(M, axis=1)
    col_norms = np.linalg.norm(M, axis=0)
    assert row_norms.shape == (3,)
    assert col_norms.shape == (4,)
    # Seedling row: sqrt(2^2 + 4^2 + 1^2 + 3^2) = sqrt(30), by hand.
    assert np.allclose(row_norms[0], np.sqrt(30.0), atol=TOL)
    # base column: sqrt(2^2 + 0^2 + 6^2) = sqrt(40), by hand.
    assert np.allclose(col_norms[0], np.sqrt(40.0), atol=TOL)


def test_meaning_three_transformation_consumes_columns_returns_rows(M):
    prices = np.array(PRICE_PER_LITRE)
    out = (M * prices).sum(axis=1)
    assert out.shape == (3,)
    assert out.tolist() == COST_PER_BAG_PENCE == [36, 27, 84]


def test_the_from_scratch_transformation_agrees(M):
    assert Matrix(RECIPES).apply_to(PRICE_PER_LITRE) == COST_PER_BAG_PENCE
    assert np.allclose(
        Matrix(RECIPES).apply_to(PRICE_PER_LITRE),
        (M * np.array(PRICE_PER_LITRE)).sum(axis=1),
        atol=TOL,
    )


def test_a_transformation_rejects_a_vector_of_the_wrong_length():
    with pytest.raises(ShapeMismatch):
        Matrix(RECIPES).apply_to([1, 2, 3])


# --------------------------------------------------------------------------
# 3. Views and copies
# --------------------------------------------------------------------------


def test_reshape_returns_a_view_and_mutation_travels(M):
    flat = M.reshape(12)
    assert flat.tolist() == [2, 4, 1, 3, 0, 5, 2, 7, 6, 1, 4, 2]
    assert np.shares_memory(M, flat)
    assert flat.base is M
    flat[0] = 99
    assert M[0, 0] == 99, "writing through the view changed the original"


def test_copy_breaks_the_link(M):
    independent = M.copy().reshape(12)
    independent[0] = 99
    assert M[0, 0] == 2
    assert not np.shares_memory(M, independent)


def test_base_is_not_the_test_for_independence(M):
    """`.base is None` is False here, and the array is still independent."""
    independent = M.copy().reshape(12)
    assert independent.base is not None, "base points at the anonymous copy"
    assert not np.shares_memory(M, independent)


def test_a_slice_is_a_view(M):
    column = M[:, 2]
    assert column.tolist() == [1, 2, 4]
    assert np.shares_memory(M, column)
    column[0] = 50
    assert M[0, 2] == 50


def test_fancy_indexing_copies(M):
    picked = M[[0, 2]]
    assert picked.shape == (2, 4)
    assert not np.shares_memory(M, picked)
    picked[0, 0] = 77
    assert M[0, 0] == 2


def test_ravel_views_and_flatten_copies(M):
    assert np.shares_memory(M, M.ravel())
    assert not np.shares_memory(M, M.flatten())
    assert M.ravel().tolist() == M.flatten().tolist()


def test_transpose_is_a_view(M):
    t = M.T
    assert np.shares_memory(M, t)
    t[0, 0] = 42
    assert M[0, 0] == 42


def test_impossible_reshape_raises(M):
    with pytest.raises(ValueError):
        M.reshape(5, 3)


def test_minus_one_infers_the_missing_dimension(M):
    assert M.reshape(6, -1).shape == (6, 2)
    assert M.reshape(-1).shape == (12,)


# --------------------------------------------------------------------------
# 4. Broadcasting
# --------------------------------------------------------------------------


def test_broadcasting_a_row_vector_across_every_row(M):
    prices = np.array(PRICE_PER_LITRE)
    scaled = M * prices
    assert scaled.shape == (3, 4)
    assert scaled[0].tolist() == [20, 8, 5, 3]
    assert scaled[1].tolist() == [0, 10, 10, 7]
    assert scaled[2].tolist() == [60, 2, 20, 2]


def test_broadcasting_copies_nothing():
    prices = np.array(PRICE_PER_LITRE)
    stretched = np.broadcast_to(prices, (3, 4))
    assert stretched.shape == (3, 4)
    assert np.shares_memory(stretched, prices)
    assert stretched.flags.writeable is False


def test_broadcasting_failure_is_a_ValueError(M):
    with pytest.raises(ValueError) as caught:
        M + np.array([100, 200, 300])
    assert "could not be broadcast together" in str(caught.value)


def test_the_failure_is_fixed_by_naming_the_axis(M):
    column = np.array([100, 200, 300]).reshape(3, 1)
    result = M + column
    assert result.shape == (3, 4)
    assert result[0].tolist() == [102, 104, 101, 103]
    assert result[2].tolist() == [306, 301, 304, 302]


def test_broadcast_shapes_predicts_both_outcomes():
    assert np.broadcast_shapes((3, 4), (4,)) == (3, 4)
    assert np.broadcast_shapes((3, 4), (3, 1)) == (3, 4)
    assert np.broadcast_shapes((3, 1), (1, 4)) == (3, 4)
    with pytest.raises(ValueError):
        np.broadcast_shapes((3, 4), (3,))


def test_the_square_matrix_trap_is_silent_and_wrong():
    S = np.array(
        [
            [1.0, 2.0, 3.0, 4.0],
            [10.0, 20.0, 30.0, 40.0],
            [100.0, 200.0, 300.0, 400.0],
            [1000.0, 2000.0, 3000.0, 4000.0],
        ]
    )
    wrong = S - S.mean(axis=1)
    right = S - S.mean(axis=1, keepdims=True)
    # Both have the same shape. Only one is what was meant.
    assert wrong.shape == right.shape == (4, 4)
    assert not np.allclose(wrong, right, atol=TOL)
    # Row-centred data must have every ROW summing to zero.
    assert np.allclose(right.sum(axis=1), np.zeros(4), atol=1e-9)
    assert not np.allclose(wrong.sum(axis=1), np.zeros(4), atol=1e-9)
    # Concretely: the silent version subtracted ROW j's mean from COLUMN j.
    # Row 0's mean is 2.5, so the whole of column 0 lost 2.5.
    assert np.allclose(S.mean(axis=1), [2.5, 25.0, 250.0, 2500.0], atol=TOL)
    assert np.allclose(wrong[2, 0], 100.0 - 2.5, atol=TOL)
    assert np.allclose(right[2, 0], 100.0 - 250.0, atol=TOL)


def test_the_same_mistake_on_a_non_square_matrix_is_loud(M):
    with pytest.raises(ValueError):
        M - M.mean(axis=1)


# --------------------------------------------------------------------------
# 5. Axes
# --------------------------------------------------------------------------


def test_axis_zero_collapses_the_rows(M):
    total = M.sum(axis=0)
    assert total.shape == (4,)
    assert total.tolist() == LITRES_PER_INGREDIENT == [8, 10, 7, 12]


def test_axis_one_collapses_the_columns(M):
    total = M.sum(axis=1)
    assert total.shape == (3,)
    assert total.tolist() == LITRES_PER_BAG == [10, 14, 13]


def test_no_axis_collapses_everything(M):
    assert M.sum() == 37
    assert M.sum(axis=0).sum() == M.sum(axis=1).sum() == 37
    assert np.shape(M.sum()) == ()


def test_means_along_each_axis(M):
    assert np.allclose(M.mean(axis=0), [8 / 3, 10 / 3, 7 / 3, 4.0], atol=TOL)
    assert np.allclose(M.mean(axis=1), [2.5, 3.5, 3.25], atol=TOL)


def test_argmax_returns_positions_not_values(M):
    assert np.argmax(M, axis=0).tolist() == [2, 1, 2, 1]
    assert np.argmax(M, axis=1).tolist() == [1, 3, 0]
    assert INGREDIENT_NAMES[np.argmax(M, axis=1)[2]] == "base"


def test_keepdims_leaves_a_one_where_the_axis_was(M):
    assert M.sum(axis=0, keepdims=True).shape == (1, 4)
    assert M.sum(axis=1, keepdims=True).shape == (3, 1)


def test_keepdims_is_what_makes_row_normalisation_work(M):
    share = M / M.sum(axis=1, keepdims=True)
    assert share.shape == (3, 4)
    assert np.allclose(share.sum(axis=1), np.ones(3), atol=TOL)
    # Seedling is 2 of 10 litres base, so exactly 0.2.
    assert np.allclose(share[0, 0], 0.2, atol=TOL)

lesson_id: D100
day: 100
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-100-matrices-and-what-they-represent
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - .venv/bin/python3 -c "import numpy; print(numpy.__version__)"
run_commands:
  - 'cd examples && ../.venv/bin/python3 01_matrix_from_scratch.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_three_meanings.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_views_and_copies.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_broadcasting.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_axes.py && cd ..'
  - .venv/bin/pytest examples -q -p no:cacheprovider
  - .venv/bin/pytest starter -q -p no:cacheprovider
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - "find . -type d -name '__pycache__' -prune -exec rm -rf -- {} +"
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 30
last_executed: '2026-08-16'
executed_on: 'macOS 26.5.2 (Apple Silicon, arm64), Python 3.14.0, numpy 2.5.2, pytest 9.1.1, bash 3.2.57 — bash tests/run_tests.sh -> 41 checks, 0 failure(s), exit 0; pytest examples -> 41 passed; pytest starter -> 1 passed, 32 skipped on an untouched checkout, and 33 passed against a fully solved copy of starter/ kept outside the lab. All five reference scripts exit 0 with every internal assertion holding. Network is needed once to install numpy and pytest; nothing else in the lab opens a socket, and section 7 of the harness greps the sources to prove it. Section 6 re-runs the harness with one expectation deliberately swapped for the wrong axis answer 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.'

# Dependencies for the Day 100 lab

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

| Package | Pinned version | Why this lab needs it |
| --- | --- | --- |
| `numpy` | `2.5.2` | The array library the lesson is about. It supplies `ndarray`, `.shape`, `.reshape`, broadcasting, the `axis` argument, and `numpy.shares_memory`, which is the tool the view-versus-copy exercise turns on. |
| `pytest` | `9.1.1` | The test runner from Days 071–074. Nothing new here except what it is pointed at. |

Nothing else is required. Both suites and all five reference scripts run on
those two packages plus the standard library.

## Why numpy is pinned and pytest almost need not be

NumPy 2.0 changed how a lone array element prints — `np.int64(36)` rather than
the bare `36` NumPy 1.x produced — and the captured files in
`../expected-output/` would not match across that boundary if the scripts
printed raw array elements. They do not: every printed array goes through
`.tolist()` first. The pin is still there because the version is *checked*
rather than assumed. Section 1 of `tests/run_tests.sh` reads the installed
version and compares it against `requirements.txt`, so a mismatch is reported
at the top of the run instead of surfacing later as a confusing diff.

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

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

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

## Licences

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

## One-time install

From the lab directory:

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

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

## Network

Installing needs the network, once. **Nothing else in this lab does.** No
script opens a socket, reads a URL or contacts a service, and section 7 of
`tests/run_tests.sh` greps every file under `examples/` and `starter/` for the
patterns that would indicate otherwise.

If you are offline and already have NumPy available somewhere, you do not need
the install at all — see the next section.

## Running without a lab-local environment

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

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

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

## What you would give up without NumPy

Exercise 1 — the from-scratch matrix class — needs nothing but the standard
library, and you can complete it on a bare `python3`. Everything after it
compares your class against NumPy, or demonstrates something (views,
broadcasting, `axis`) that only exists because NumPy exists. Those parts
cannot be faked, and the lab does not pretend otherwise.

numpy==2.5.2
pytest==9.1.1

# The brief — The Same Numbers, Three Ways

Work through the five exercises below in order. Everything you write goes into
two files in this directory: `matrix.py` (exercise 1) and `answers.py`
(exercises 2 to 5). Nothing else needs editing.

Run this after every change, from the **lab directory**, one level up:

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

Anything you have not attempted is **skipped**, not failed, so the summary line
is a running score. On an untouched checkout it reads `1 passed, 32 skipped`.
When it reads `33 passed`, you are finished.

## The data

Three invented potting mixes from an invented garden centre, described by four
ingredients, in litres per bag:

|  | base | bark | grit | compost |
| --- | --- | --- | --- | --- |
| **Seedling** | 2 | 4 | 1 | 3 |
| **Container** | 0 | 5 | 2 | 7 |
| **Alpine** | 6 | 1 | 4 | 2 |

and the ingredient prices, in pence per litre: base 10, bark 2, grit 5,
compost 1.

Twelve numbers, chosen so that every answer in this lab can be worked out on
paper in under a minute. That is the only property that matters here. The
garden centre, the recipes and the prices are all invented.

Three rows and four columns, deliberately: because 3 and 4 are different
numbers, the **shape of an answer tells you which operation produced it**. On a
square matrix it would not, and exercise 4 shows you what that costs.

---

## Exercise 1 — build the matrix yourself (`matrix.py`)

Fill in the six methods marked `EXERCISE` in `matrix.py`:

| Method | What it must do |
| --- | --- |
| `shape` | Return `(rows, columns)` — rows first |
| `__getitem__` | Support `m[i, j]` from 0; `TypeError` for a non-pair, `IndexError` naming the shape for an out-of-range or negative index |
| `transpose` | Swap rows and columns: `(r, c)` becomes `(c, r)` |
| `add` | Elementwise, identical shapes only; `ShapeMismatch` for a size clash, `TypeError` for a plain list |
| `scale` | Multiply every entry by one number |
| `identity` | The `n` by `n` matrix with 1 on the diagonal and 0 elsewhere |

Do not import numpy in that file. The whole value of the exercise is that
nothing is done for you; the tests then check your class against NumPy, which
is the point at which you find out whether you agreed with it.

Ten of the thirty-three tests belong to this exercise.

## Exercise 2 — the three meanings (`answers.py`, section 2)

Six predictions. Read the same twelve numbers as a **table** (which row is
Alpine, which column is grit), then as a **transformation** (apply the matrix
to the price vector and you get the cost of one bag of each mix).

Work the three costs out on paper. Each is four multiplications and three
additions. Then note which length the answer came out as, and which length was
consumed.

## Exercise 3 — views and copies (`answers.py`, section 3)

Five predictions about whether writing through a second name changes the
first. Answer them from what you believe *before* running anything — this is
the section where most people discover their mental model was wrong, and that
discovery only happens if you commit to an answer.

## Exercise 4 — broadcasting (`answers.py`, section 4)

Five predictions. Apply the rule by hand, right to left, before you let NumPy
apply it for you:

1. Line the two shapes up from the **right-hand** end.
2. A missing entry on the left of the shorter shape counts as 1.
3. Two dimensions are compatible when they are equal, or one of them is 1.
4. If any pair is neither, the operation is an error.
5. The result takes the larger of each pair.

## Exercise 5 — axis=0 against axis=1 (`answers.py`, section 5)

Seven predictions, and the one rule that settles it permanently:

> **The axis you name is the axis that disappears.**

A `(3, 4)` array summed with `axis=0` loses the 3 and returns shape `(4,)`. The
same array summed with `axis=1` loses the 4 and returns shape `(3,)`.

---

## When you are finished

Compare your `matrix.py` with `../examples/matrix.py`. They should agree on
behaviour; they need not agree on wording. Then read
`../examples/03_views_and_copies.py` and `../examples/04_broadcasting.py`,
which go several steps past what the predictions asked for.

"""Exercises 2 to 5 — predict first, then let NumPy tell you.

Every value below is a prediction you make BEFORE running anything. Replace
each `None` with your answer, then check them from the lab directory:

    .venv/bin/pytest starter -q

Each prediction still set to None is skipped rather than failed, so the test
run is a running score. A wrong prediction fails with both numbers printed.

Predicting first matters more than it sounds. Broadcasting and axis arguments
are where NumPy stops doing the obvious thing, and the only way to find out
whether your mental model is right is to commit to an answer while it can
still be wrong.

The matrix everything below refers to — three invented potting mixes described
by four ingredients, in litres per bag:

                base   bark   grit   compost
    Seedling       2      4      1         3
    Container      0      5      2         7
    Alpine         6      1      4         2

and the ingredient prices, in pence per litre: base 10, bark 2, grit 5,
compost 1.
"""

# ---------------------------------------------------------------------------
# Exercise 2 — shape and the three meanings
# ---------------------------------------------------------------------------

# 2.1 The shape of the matrix above, as a (rows, columns) tuple.
SHAPE_OF_M = None

# 2.2 The shape of its transpose.
SHAPE_OF_M_T = None

# 2.3 Read as a TABLE: the Alpine row, as a plain list of four integers.
ALPINE_ROW = None

# 2.4 Read as a TABLE: the grit column, as a plain list of three integers.
GRIT_COLUMN = None

# 2.5 Read as a TRANSFORMATION: apply the matrix to the price vector and you
#     get the cost of one bag of each mix, in pence. Work all three out on
#     paper — each is four multiplications and three additions — and write
#     them here as a list of three integers, in the mix order above.
COST_PER_BAG_PENCE = None

# 2.6 A (3, 4) matrix applied to a vector returns a vector of what length?
LENGTH_OF_TRANSFORMED_VECTOR = None

# ---------------------------------------------------------------------------
# Exercise 3 — views and copies
# ---------------------------------------------------------------------------

# 3.1 M is the (3, 4) array. You write:
#         flat = M.reshape(12)
#         flat[0] = 99
#     What is M[0, 0] afterwards? An integer.
M_00_AFTER_WRITING_THROUGH_RESHAPE = None

# 3.2 Same again, but:
#         independent = M.copy().reshape(12)
#         independent[0] = 99
#     What is M[0, 0] afterwards? An integer.
M_00_AFTER_WRITING_THROUGH_A_COPY = None

# 3.3 Does a basic slice — M[:, 2] — share memory with M? True or False.
SLICE_SHARES_MEMORY = None

# 3.4 Does fancy indexing — M[[0, 2]] — share memory with M? True or False.
FANCY_INDEX_SHARES_MEMORY = None

# 3.5 Does M.T share memory with M? True or False.
TRANSPOSE_SHARES_MEMORY = None

# ---------------------------------------------------------------------------
# Exercise 4 — broadcasting
# ---------------------------------------------------------------------------

# 4.1 Apply the rule by hand: shape (3, 4) combined with shape (4,).
#     Write the resulting shape as a tuple, or the string "error" if the rule
#     rejects it.
BROADCAST_3x4_WITH_4 = None

# 4.2 Shape (3, 4) combined with shape (3,). Tuple, or "error".
BROADCAST_3x4_WITH_3 = None

# 4.3 Shape (3, 4) combined with shape (3, 1). Tuple, or "error".
BROADCAST_3x4_WITH_3x1 = None

# 4.4 Shape (3, 1) combined with shape (1, 4). Tuple, or "error".
BROADCAST_3x1_WITH_1x4 = None

# 4.5 When broadcasting is rejected, which exception class does NumPy raise?
#     Write the class itself, not its name as a string — for example
#     BROADCAST_FAILURE_EXCEPTION = KeyError
BROADCAST_FAILURE_EXCEPTION = None

# ---------------------------------------------------------------------------
# Exercise 5 — axis=0 against axis=1
# ---------------------------------------------------------------------------

# 5.1 M.sum(axis=0) — write the resulting shape as a tuple.
SHAPE_OF_SUM_AXIS_0 = None

# 5.2 M.sum(axis=1) — write the resulting shape as a tuple.
SHAPE_OF_SUM_AXIS_1 = None

# 5.3 "How many litres are in each bag?" is one number per mix. Which axis
#     argument answers it — 0 or 1?
AXIS_FOR_LITRES_PER_BAG = None

# 5.4 "How many litres of each ingredient does one bag of every mix need?" is
#     one number per ingredient. Which axis argument answers it — 0 or 1?
AXIS_FOR_LITRES_PER_INGREDIENT = None

# 5.5 Work out both totals on paper and write them as lists of integers.
LITRES_PER_BAG = None
LITRES_PER_INGREDIENT = None

# 5.6 M.sum(axis=1, keepdims=True) — the resulting shape, as a tuple.
SHAPE_OF_SUM_AXIS_1_KEEPDIMS = None

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

Both `examples/` and `starter/` contain a module called `matrix`, 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
`matrix` 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
`matrix` or `dataset` 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 ("matrix", "dataset"):
    module = sys.modules.get(name)
    origin = getattr(module, "__file__", "") or ""
    if module is not None and not origin.startswith(HERE):
        del sys.modules[name]

"""Exercise 1 — build the matrix yourself, on nothing but nested lists.

Fill in the five methods marked EXERCISE below. Each one is a few lines. The
point is not difficulty; it is that after writing them you will never again be
unsure what `.shape`, `.T` or an elementwise sum is actually doing, because you
will have done it.

Check your work as you go, from the LAB DIRECTORY (one level up from here):

    .venv/bin/pytest starter -q

Every test for a method you have not written yet is skipped, not failed, so the
output is a running score rather than a wall of red. When all of them pass,
compare your file with examples/matrix.py — they should agree on behaviour,
not necessarily on wording.

Rules for all five: return a NEW Matrix, never modify self, and do not import
numpy in this file. The whole value of the exercise is that nothing here is
done for you.
"""

from __future__ import annotations


class ShapeMismatch(ValueError):
    """Raised when two matrices cannot be combined because their shapes differ.

    Subclassing ValueError is deliberate: NumPy raises ValueError for the same
    situation, so one `except ValueError` catches your class and NumPy alike.
    """


class Matrix:
    """A rectangular grid of numbers, stored as a list of row lists."""

    def __init__(self, rows):
        rows = [list(row) for row in rows]
        if not rows or not rows[0]:
            raise ValueError("a matrix needs at least one row and one column")
        width = len(rows[0])
        for i, row in enumerate(rows):
            if len(row) != width:
                raise ValueError(
                    f"row {i} has {len(row)} entries but row 0 has {width}; "
                    "a matrix is rectangular"
                )
        self._rows = rows

    # ---------------------------------------------------------------- 1.1 --
    @property
    def shape(self):
        """EXERCISE 1.1 — return (number of rows, number of columns).

        Rows first. That order is a convention rather than a law, and it is
        worth saying out loud once so you never have to wonder again.

        Hint: len(self._rows) and len(self._rows[0]).
        """
        raise NotImplementedError("write Matrix.shape")

    # ---------------------------------------------------------------- 1.2 --
    def __getitem__(self, position):
        """EXERCISE 1.2 — support m[i, j], counting rows and columns from 0.

        `position` arrives as a tuple (i, j). Two requirements beyond the
        obvious lookup:

          * if `position` is not a 2-tuple, raise TypeError with a message
            that says how to index a Matrix;
          * if i or j is outside the matrix, raise IndexError with a message
            that includes the string f"shape {self.shape}" — the tests check
            for exactly that, because an error that does not tell you the
            shape makes you go and print it yourself.

        Note that negative indices should be rejected too: Python lists accept
        m[-1], and a matrix that silently wraps around is a bug generator.
        """
        raise NotImplementedError("write Matrix.__getitem__")

    def row(self, i):
        """Row i as a plain list. Already written, and it copies on purpose."""
        return list(self._rows[i])

    def col(self, j):
        """Column j as a plain list. Note that this has to walk every row."""
        return [row[j] for row in self._rows]

    # ---------------------------------------------------------------- 1.3 --
    def transpose(self):
        """EXERCISE 1.3 — swap rows and columns; (r, c) becomes (c, r).

        Entry (i, j) of the result is entry (j, i) of the original. A nested
        list comprehension does it in one line, but write it as two loops
        first if that reads more clearly to you.
        """
        raise NotImplementedError("write Matrix.transpose")

    @property
    def T(self):
        """Spelled the way NumPy spells it. Already written."""
        return self.transpose()

    # ---------------------------------------------------------------- 1.4 --
    def add(self, other):
        """EXERCISE 1.4 — elementwise addition, with no broadcasting at all.

          * if `other` is not a Matrix, raise TypeError, and say in the message
            that this class has no broadcasting;
          * if the shapes differ, raise ShapeMismatch naming both shapes;
          * otherwise return a new Matrix of the entrywise sums.
        """
        raise NotImplementedError("write Matrix.add")

    def __add__(self, other):
        return self.add(other)

    # ---------------------------------------------------------------- 1.5 --
    def scale(self, k):
        """EXERCISE 1.5 — multiply every entry by the single number k."""
        raise NotImplementedError("write Matrix.scale")

    def __mul__(self, k):
        return self.scale(k)

    __rmul__ = __mul__

    # ---------------------------------------------------------------- 1.6 --
    @classmethod
    def identity(cls, n):
        """EXERCISE 1.6 — the n by n matrix with 1 on the diagonal, 0 elsewhere.

        This is the matrix that leaves every vector exactly as it found it, and
        recognising it on sight is worth more than it looks.
        """
        raise NotImplementedError("write Matrix.identity")

    # -- already written, so you have something to test against -------------

    @classmethod
    def zeros(cls, n_rows, n_cols):
        return cls([[0] * n_cols for _ in range(n_rows)])

    def apply_to(self, vector):
        """Treat the matrix as a transformation and apply it to one vector.

        Each output entry is one row multiplied entry by entry against the
        input vector and then summed. An (r, c) matrix therefore eats a vector
        of length c and returns a vector of length r.

        Written for you because the packed name for it is Day 101's subject.
        """
        vector = list(vector)
        n_rows, n_cols = self.shape
        if len(vector) != n_cols:
            raise ShapeMismatch(
                f"a {self.shape} matrix transforms a vector of length {n_cols}, "
                f"but this vector has length {len(vector)}"
            )
        return [
            sum(self[i, j] * vector[j] for j in range(n_cols)) for i in range(n_rows)
        ]

    def to_lists(self):
        return [list(row) for row in self._rows]

    def __eq__(self, other):
        if not isinstance(other, Matrix):
            return NotImplemented
        return self._rows == other._rows

    def __repr__(self):
        body = ", ".join(repr(row) for row in self._rows)
        return f"Matrix([{body}])"

    def format(self, width=4):
        return "\n".join(
            " ".join(f"{value:>{width}}" for value in row) for row in self._rows
        )

"""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 numpy.allclose with atol=TOL, stated below.
"""

import numpy as np
import pytest

import answers
from matrix import Matrix, ShapeMismatch

TOL = 1e-12

RECIPES = [
    [2, 4, 1, 3],
    [0, 5, 2, 7],
    [6, 1, 4, 2],
]
PRICE_PER_LITRE = [10, 2, 5, 1]


def written(fn, *args, **kwargs):
    """Run part of your Matrix, 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


@pytest.fixture
def M():
    return np.array(RECIPES)


@pytest.fixture
def A():
    return Matrix(RECIPES)


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


def test_0_the_environment_is_ready():
    """Always passes once the install worked. Everything below is your work."""
    assert np.__version__, "numpy is importable"
    # The constructor is written for you, and it already refuses a ragged grid:
    # a matrix is rectangular, and that is not a detail you get to skip.
    with pytest.raises(ValueError):
        Matrix([[1, 2, 3], [4, 5]])


# -- Exercise 1: your Matrix ------------------------------------------------


def test_1_1_shape(A):
    assert written(lambda: A.shape) == (3, 4)


def test_1_2_indexing(A):
    assert written(lambda: A[0, 0]) == 2
    assert A[1, 3] == 7
    assert A[2, 3] == 2


def test_1_2_index_errors_name_the_shape(A):
    written(lambda: A[0, 0])
    with pytest.raises(IndexError) as caught:
        A[3, 0]
    assert "shape (3, 4)" in str(caught.value)
    with pytest.raises(IndexError):
        A[0, -1]


def test_1_2_indexing_with_one_number_is_a_TypeError(A):
    written(lambda: A[0, 0])
    with pytest.raises(TypeError):
        A[0]


def test_1_3_transpose(A, M):
    assert written(lambda: A.T.shape) == (4, 3)
    assert A.T.to_lists() == M.T.tolist()
    assert A.T.T == A


def test_1_4_addition(A, M):
    B = Matrix([[1, 0, 0, 1], [2, 2, 2, 2], [0, 3, 0, 3]])
    assert written(lambda: (A + B).to_lists()) == (M + np.array(B.to_lists())).tolist()


def test_1_4_addition_rejects_a_shape_mismatch(A):
    written(lambda: A.add(A))
    with pytest.raises(ShapeMismatch):
        A.add(Matrix([[1, 2], [3, 4]]))
    with pytest.raises(ValueError):
        A.add(Matrix([[1, 2], [3, 4]]))


def test_1_4_addition_refuses_a_plain_list(A):
    written(lambda: A.add(A))
    with pytest.raises(TypeError) as caught:
        A.add([100, 200, 300, 400])
    assert "broadcasting" in str(caught.value)


def test_1_5_scale(A, M):
    assert written(lambda: (A * 3).to_lists()) == (M * 3).tolist()
    assert (3 * A).to_lists() == (M * 3).tolist()


def test_1_6_identity():
    assert written(Matrix.identity, 3).to_lists() == np.eye(3, dtype=int).tolist()
    v = [1.5, -2.0, 0.25]
    assert np.allclose(Matrix.identity(3).apply_to(v), v, atol=TOL)


# -- Exercise 2: shape and the three meanings -------------------------------


def test_2_1_shape_of_M(M):
    assert predicted("SHAPE_OF_M") == M.shape


def test_2_2_shape_of_transpose(M):
    assert predicted("SHAPE_OF_M_T") == M.T.shape


def test_2_3_alpine_row(M):
    assert predicted("ALPINE_ROW") == M[2].tolist()


def test_2_4_grit_column(M):
    assert predicted("GRIT_COLUMN") == M[:, 2].tolist()


def test_2_5_cost_per_bag(M):
    real = (M * np.array(PRICE_PER_LITRE)).sum(axis=1)
    assert predicted("COST_PER_BAG_PENCE") == real.tolist()


def test_2_6_transformed_vector_length(M):
    real = (M * np.array(PRICE_PER_LITRE)).sum(axis=1)
    assert predicted("LENGTH_OF_TRANSFORMED_VECTOR") == len(real)


# -- Exercise 3: views and copies -------------------------------------------


def test_3_1_writing_through_a_reshape(M):
    guess = predicted("M_00_AFTER_WRITING_THROUGH_RESHAPE")
    flat = M.reshape(12)
    flat[0] = 99
    assert guess == M[0, 0]


def test_3_2_writing_through_a_copy(M):
    guess = predicted("M_00_AFTER_WRITING_THROUGH_A_COPY")
    independent = M.copy().reshape(12)
    independent[0] = 99
    assert guess == M[0, 0]


def test_3_3_slice_shares_memory(M):
    assert predicted("SLICE_SHARES_MEMORY") == np.shares_memory(M, M[:, 2])


def test_3_4_fancy_index_shares_memory(M):
    assert predicted("FANCY_INDEX_SHARES_MEMORY") == np.shares_memory(M, M[[0, 2]])


def test_3_5_transpose_shares_memory(M):
    assert predicted("TRANSPOSE_SHARES_MEMORY") == np.shares_memory(M, M.T)


# -- Exercise 4: broadcasting -----------------------------------------------


def _outcome(left, right):
    try:
        return np.broadcast_shapes(left, right)
    except ValueError:
        return "error"


@pytest.mark.parametrize(
    "name, left, right",
    [
        ("BROADCAST_3x4_WITH_4", (3, 4), (4,)),
        ("BROADCAST_3x4_WITH_3", (3, 4), (3,)),
        ("BROADCAST_3x4_WITH_3x1", (3, 4), (3, 1)),
        ("BROADCAST_3x1_WITH_1x4", (3, 1), (1, 4)),
    ],
)
def test_4_broadcast_outcomes(name, left, right):
    assert predicted(name) == _outcome(left, right)


def test_4_5_failure_exception_class(M):
    guess = predicted("BROADCAST_FAILURE_EXCEPTION")
    with pytest.raises(guess):
        M + np.array([100, 200, 300])


# -- Exercise 5: axes -------------------------------------------------------


def test_5_1_shape_of_sum_axis_0(M):
    assert predicted("SHAPE_OF_SUM_AXIS_0") == M.sum(axis=0).shape


def test_5_2_shape_of_sum_axis_1(M):
    assert predicted("SHAPE_OF_SUM_AXIS_1") == M.sum(axis=1).shape


def test_5_3_axis_for_litres_per_bag(M):
    axis = predicted("AXIS_FOR_LITRES_PER_BAG")
    assert M.sum(axis=axis).tolist() == [10, 14, 13]


def test_5_4_axis_for_litres_per_ingredient(M):
    axis = predicted("AXIS_FOR_LITRES_PER_INGREDIENT")
    assert M.sum(axis=axis).tolist() == [8, 10, 7, 12]


def test_5_5_the_two_totals(M):
    assert predicted("LITRES_PER_BAG") == M.sum(axis=1).tolist()
    assert predicted("LITRES_PER_INGREDIENT") == M.sum(axis=0).tolist()


def test_5_6_keepdims_shape(M):
    assert predicted("SHAPE_OF_SUM_AXIS_1_KEEPDIMS") == M.sum(axis=1, keepdims=True).shape

#!/usr/bin/env bash
# Tests for the Day 100 lab. Run from the lab directory:
#   bash tests/run_tests.sh
#
# The harness proves six specific claims the lesson makes, and it proves each
# one by running code and reading real values rather than by reading source:
#
#   * the from-scratch matrix agrees with NumPy on shape, indexing, transpose,
#     addition and scalar multiplication — and refuses the one thing NumPy
#     does for free, which is broadcasting;
#   * the same twelve numbers answer three different questions according to
#     whether they are read as a table, as vectors, or as a transformation;
#   * a reshape is a VIEW: writing through it changes the original, and
#     .copy() breaks that link;
#   * broadcasting succeeds for (3, 4) with (4,) and fails for (3, 4) with
#     (3,) with a ValueError naming both shapes — the exact exception type is
#     asserted, not merely "it raised something";
#   * axis=0 collapses the rows and axis=1 collapses the columns, checked
#     against totals small enough to work out on paper;
#   * nothing is left behind on disk.
#
# Everything runs offline. Nothing binds a port, nothing writes outside the
# lab, nothing needs a key. Deterministic, non-interactive, exits 0 only if
# every check passes.
set -u

export PYTHONDONTWRITEBYTECODE=1

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

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

failures=0
checks=0

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

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

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

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

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

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

echo "Day 100 — The Same Numbers, Three Ways"
echo

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

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

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

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

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

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

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

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

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 40 ]; then
  check "the reference suite ran at least 40 tests (ran ${ref_passed})" "yes"
else
  check "the reference suite ran at least 40 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 a module called `matrix`, and
# pytest imports test files by putting their directory on sys.path — so
# collecting both suites at once would otherwise let the starter tests import
# the REFERENCE solution and report unwritten exercises as passing. Each
# directory's conftest.py prevents that. This check proves it still does:
# across both suites, the skip count must be unchanged.
both_out="$(cd "${lab_dir}" && "${pytest_bin}" -q -p no:cacheprovider 2>&1)"
start_skipped="$(printf '%s\n' "${start_out}" | grep -o '[0-9][0-9]* skipped' | head -1 | cut -d' ' -f1)"
both_skipped="$(printf '%s\n' "${both_out}" | grep -o '[0-9][0-9]* skipped' | head -1 | cut -d' ' -f1)"
check_eq "collecting both suites at once does not turn skips into passes" \
  "${start_skipped:-none}" "${both_skipped:-none}"

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

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

from dataset import PRICE_PER_LITRE, RECIPES
from matrix import Matrix

M = np.array(RECIPES)
A = Matrix(RECIPES)

print("shape", M.shape, A.shape)
print("transpose_shape", M.T.shape)
print("cost", (M * np.array(PRICE_PER_LITRE)).sum(axis=1).tolist())
print("cost_scratch", A.apply_to(PRICE_PER_LITRE))
print("axis0", M.sum(axis=0).tolist(), M.sum(axis=0).shape)
print("axis1", M.sum(axis=1).tolist(), M.sum(axis=1).shape)
print("keepdims", M.sum(axis=1, keepdims=True).shape)

flat = M.reshape(12)
flat[0] = 99
print("reshape_is_view", M[0, 0], np.shares_memory(M, flat))

N = np.array(RECIPES)
independent = N.copy().reshape(12)
independent[0] = 99
print("copy_breaks_link", N[0, 0], np.shares_memory(N, independent))

P = np.array(RECIPES)
column = P[:, 2]
column[0] = 50
print("slice_is_view", P[0, 2], np.shares_memory(P, column))
print("fancy_is_copy", np.shares_memory(P, P[[0, 2]]))

print("broadcast_ok", (np.array(RECIPES) * np.array(PRICE_PER_LITRE))[0].tolist())
try:
    np.array(RECIPES) + np.array([100, 200, 300])
except Exception as exc:  # deliberately broad: the TYPE is what is asserted
    print("broadcast_fail", type(exc).__name__, "could not be broadcast together" in str(exc))
else:
    print("broadcast_fail", "NOTHING_RAISED", False)

# The silent trap, on a square matrix, where no exception is raised at all.
S = np.array([[1.0, 2, 3, 4], [10, 20, 30, 40], [100, 200, 300, 400], [1000, 2000, 3000, 4000]])
wrong = S - S.mean(axis=1)
right = S - S.mean(axis=1, keepdims=True)
print("trap_same_shape", wrong.shape == right.shape)
print("trap_row_sums_zero", bool(np.allclose(right.sum(axis=1), 0, atol=1e-9)),
      bool(np.allclose(wrong.sum(axis=1), 0, atol=1e-9)))

try:
    A.add([100, 200, 300, 400])
except TypeError:
    print("scratch_refuses_broadcast", "TypeError")
else:
    print("scratch_refuses_broadcast", "NOTHING_RAISED")
PY
)"

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

check_eq "shape is (3, 4) in both implementations" "(3, 4) (3, 4)" "$(get shape)"
check_eq "transpose is (4, 3)" "(4, 3)" "$(get transpose_shape)"
check_eq "the transformation returns three costs in pence" "[36, 27, 84]" "$(get cost)"
check_eq "the from-scratch transformation returns the same three" "[36, 27, 84]" "$(get cost_scratch)"
# Section 6 re-runs this script with D100_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_axis0="[8, 10, 7, 12] (4,)"
if [ -n "${D100_SELF_TEST:-}" ]; then
  expected_axis0="[10, 14, 13] (3,)"   # the axis=1 answer, deliberately wrong here
fi
check_eq "axis=0 collapses the rows and returns four numbers" "${expected_axis0}" "$(get axis0)"
check_eq "axis=1 collapses the columns and returns three numbers" "[10, 14, 13] (3,)" "$(get axis1)"
check_eq "keepdims leaves a 1 where the axis was" "(3, 1)" "$(get keepdims)"
check_eq "writing through a reshape changed the original" "99 True" "$(get reshape_is_view)"
check_eq "writing through a copy did not" "2 False" "$(get copy_breaks_link)"
check_eq "writing through a slice changed the original" "50 True" "$(get slice_is_view)"
check_eq "fancy indexing gave a copy, not a view" "False" "$(get fancy_is_copy)"
check_eq "(3, 4) times (4,) broadcasts across every row" "[20, 8, 5, 3]" "$(get broadcast_ok)"
check_eq "(3, 4) plus (3,) raises ValueError naming both shapes" "ValueError True" "$(get broadcast_fail)"
check_eq "the square-matrix trap produces the SAME shape as the correct answer" "True" "$(get trap_same_shape)"
check_eq "only the keepdims version leaves every row summing to zero" "True False" "$(get trap_row_sums_zero)"
check_eq "the from-scratch class refuses to broadcast" "TypeError" "$(get scratch_refuses_broadcast)"

# --------------------------------------------------------------------------
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 WRONG axis answer, and asserts that the re-run reports the failure
# and exits non-zero. If this section passes, section 5 is not decorative.
if [ -z "${D100_SELF_TEST:-}" ]; then
  self_out="$(D100_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: axis=0 collapses the rows"*)
      check "the failing check is named in the output with both values" "yes" ;;
    *) check "the failing check is named in the output with both values" "no" ;;
  esac
  case "${self_out}" in
    *", 1 failure(s)."*)
      check "the summary line counts exactly one failure" "yes" ;;
    *) check "the summary line counts exactly one failure" "no" ;;
  esac
else
  echo "  (self-test run: section 6 does not recurse)"
fi

# --------------------------------------------------------------------------
echo
echo "7. Nothing was left behind"
# --------------------------------------------------------------------------

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

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

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

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