Hands-on lab — Day 107: Norms, Distances, and Similarity Measures

"""One query, three candidates, three measures -- and three different winners.

Run from the `examples/` directory:

    ../.venv/bin/python3 01_three_measures_three_winners.py

This is the disagreement the whole day is built on. None of the three answers
is a bug and none of them is wrong. They are answers to different questions,
and until you have decided which question you are asking, the ranking your
retrieval system returns is decided by whichever measure someone left in place.
"""

from __future__ import annotations

import catalogue
import measures
from measures import TOL

Q = catalogue.QUERY
A = catalogue.ARTICLES


def table(title: str, rows: list[tuple[str, str]]) -> None:
    print(f"  {title}")
    for left, right in rows:
        print(f"    {left:<14} {right}")
    print()


print("=" * 74)
print("1. The data")
print("=" * 74)
print()
print("  A help-centre search. Four terms, counted. The query is the reader's")
print("  own note; the three candidates are articles.")
print()
print(f"    {'':<12}" + "".join(f"{t:>10}" for t in catalogue.TERMS)
      + f"{'total':>10}")
print(f"    {'query':<12}" + "".join(f"{c:>10}" for c in Q)
      + f"{sum(Q):>10}")
for name, vec in A.items():
    print(f"    {name:<12}" + "".join(f"{c:>10}" for c in vec)
          + f"{sum(vec):>10}")
print()
print("  Cartogram is the query's exact profile at three times the length:")
print(f"    query * 3 = {tuple(3 * c for c in Q)} = {A['Cartogram']}")
assert tuple(3 * c for c in Q) == A["Cartogram"]
print()

print("=" * 74)
print("2. The three measures, computed")
print("=" * 74)
print()

l1 = {n: measures.l1_distance(Q, v) for n, v in A.items()}
l2 = {n: measures.l2_distance(Q, v) for n, v in A.items()}
linf = {n: measures.linf_distance(Q, v) for n, v in A.items()}
cos = {n: measures.cosine_similarity(Q, v) for n, v in A.items()}

header = f"    {'':<12}{'L1':>12}{'L2':>12}{'L-inf':>12}{'cosine':>12}"
print(header)
print("    " + "-" * (len(header) - 4))
for name in A:
    print(f"    {name:<12}{l1[name]:>12.4f}{l2[name]:>12.4f}"
          f"{linf[name]:>12.4f}{cos[name]:>12.4f}")
print()
print("  L1 and L2 and L-infinity are DISTANCES: smaller is better.")
print("  Cosine is a SIMILARITY: larger is better. Mixing the two up returns")
print("  the worst match with complete confidence and no error message.")
print()

print("=" * 74)
print("3. The winners")
print("=" * 74)
print()

winners = {
    "L1 (Manhattan)": measures.winner(Q, A, measures.l1_distance),
    "L2 (Euclidean)": measures.winner(Q, A, measures.l2_distance),
    "L-inf (Chebyshev)": measures.winner(Q, A, measures.linf_distance),
    "cosine similarity": measures.winner(Q, A, measures.cosine_similarity,
                                         higher_is_better=True),
}
for measure_name, best in winners.items():
    print(f"    {measure_name:<20} picks  {best}")
print()

assert winners["L1 (Manhattan)"] == "Aisle"
assert winners["L2 (Euclidean)"] == "Beacon"
assert winners["L-inf (Chebyshev)"] == "Beacon"
assert winners["cosine similarity"] == "Cartogram"
assert len({winners["L1 (Manhattan)"], winners["L2 (Euclidean)"],
            winners["cosine similarity"]}) == 3

print("  Three of the four measures name three different articles, on the")
print("  same four numbers, with no randomness and nothing to tune.")
print()

print("=" * 74)
print("4. Why each one is right, on its own terms")
print("=" * 74)
print()

table("L1 picks Aisle -- total disagreement is what it measures.", [
    ("Aisle", f"|4-4| + |3-3| + |2-2| + |1-6| = {l1['Aisle']:.0f}"),
    ("Beacon", f"|4-6| + |3-1| + |2-4| + |1-1| = {l1['Beacon']:.0f}"),
    ("", "5 is less than 6. Aisle wins, and it is a fair answer:"),
    ("", "there is genuinely less disagreement in total."),
])

table("L2 picks Beacon -- squaring punishes one BIG disagreement.", [
    ("Aisle", f"sqrt(0 + 0 + 0 + 25) = {l2['Aisle']:.4f}"),
    ("Beacon", f"sqrt(4 + 4 + 4 + 0)  = {l2['Beacon']:.4f}"),
    ("", "Aisle's single error of 5 costs 25. Beacon's three errors"),
    ("", "of 2 cost 4 each, 12 in total, even though they add to 6."),
])

table("L-infinity agrees with L2 here, for a different reason.", [
    ("Aisle", f"max(0, 0, 0, 5) = {linf['Aisle']:.0f}"),
    ("Beacon", f"max(2, 2, 2, 0) = {linf['Beacon']:.0f}"),
    ("", "It looks only at the worst single term and ignores the rest"),
    ("", "entirely. Beacon's worst is 2; Aisle's worst is 5."),
])

table("Cosine picks Cartogram -- it is scored on direction alone.", [
    ("Cartogram", f"cosine = {cos['Cartogram']:.6f}"),
    ("", "Exactly 1.0: same mix of terms, three times the length."),
    ("", f"Its L2 distance is {l2['Cartogram']:.4f}, the worst here by far,"),
    ("", "because length is the one thing cosine throws away."),
])

assert abs(cos["Cartogram"] - 1.0) <= TOL
assert l2["Cartogram"] == max(l2.values())
assert l1["Cartogram"] == max(l1.values())

print("=" * 74)
print("5. The same numbers, from NumPy")
print("=" * 74)
print()
print("  Nothing in measures.py uses NumPy. numpy.linalg.norm with an `ord`")
print("  argument IS the p-norm family, so it is an independent answer rather")
print("  than a restatement.")
print()

import numpy as np  # noqa: E402  (imported here to make the point above)

worst = 0.0
print(f"    {'':<12}{'ord=1':>12}{'ord=2':>12}{'ord=inf':>12}")
for name, vec in A.items():
    d = np.asarray(Q, dtype=float) - np.asarray(vec, dtype=float)
    row = [float(np.linalg.norm(d, ord=o)) for o in (1, 2, np.inf)]
    mine = [l1[name], l2[name], linf[name]]
    worst = max(worst, max(abs(x - y) for x, y in zip(row, mine)))
    print(f"    {name:<12}" + "".join(f"{v:>12.4f}" for v in row))
print()
print(f"  worst disagreement with measures.py: {worst:.3e}")
print(f"  stated tolerance:                    {TOL:.3e}")
assert worst <= TOL
print()

print("=" * 74)
print("6. The point")
print("=" * 74)
print()
print("  If this were a search box, the article at the top of the results is")
print("  decided by a single argument that nobody in the room has discussed.")
print()
print("  The right answer depends entirely on what you meant:")
print("    * 'closest counts overall'      -> L1  -> Aisle")
print("    * 'no single term badly wrong'  -> L2  -> Beacon")
print("    * 'same topic, any length'      -> cos -> Cartogram")
print()
print("  For text retrieval the third is almost always what you meant, which")
print("  is why cosine is the default in every vector database. The lesson to")
print("  keep is not 'use cosine'. It is that the default is a decision, and")
print("  somebody has to make it on purpose.")
print()
print("01_three_measures_three_winners.py: every assertion held.")

"""The p-norm family, and the picture that makes it click: three unit balls.

Run from the `examples/` directory:

    ../.venv/bin/python3 02_the_p_norm_family.py

L1, L2 and L-infinity are not three unrelated ideas. They are one formula with
one dial turned to three settings, and the clearest way to see the difference
is to draw the set of points each one calls "distance 1 from the origin".
A diamond, a circle, a square.
"""

from __future__ import annotations

import math

import numpy as np

import catalogue
import measures
from measures import TOL

V = (3.0, 4.0)

print("=" * 74)
print("1. One formula, one dial")
print("=" * 74)
print()
print("    p_norm(v, p) = (sum of |x| ** p) ** (1 / p)")
print()
print(f"  On the vector v = {V}, which is the 3-4-5 triangle, so the")
print("  p = 2 answer is a whole number and you can check it in your head:")
print()
print(f"    {'p':>8}{'||v||_p':>14}   note")
print("    " + "-" * 52)
notes = {
    1: "sum of |x|: 3 + 4",
    2: "the one geometry gives you",
    math.inf: "the largest single component",
}
sweep = [1, 1.5, 2, 3, 4, 8, 16, 64, math.inf]
values = []
for p in sweep:
    n = measures.p_norm(V, p)
    values.append(n)
    label = "inf" if math.isinf(p) else f"{p:g}"
    print(f"    {label:>8}{n:>14.6f}   {notes.get(p, '')}".rstrip())
print()

assert measures.p_norm(V, 1) == 7.0
assert measures.p_norm(V, 2) == 5.0
assert measures.p_norm(V, math.inf) == 4.0

print("  The value falls as p rises, and never below the largest component.")
print("  That is not a coincidence: raising a bigger number to a bigger power")
print("  makes it dominate the sum, so in the limit only the biggest survives.")
print()
for earlier, later in zip(values, values[1:]):
    assert later <= earlier + TOL, (earlier, later)
assert all(v >= measures.linf_norm(V) - TOL for v in values)
assert abs(values[-1] - measures.linf_norm(V)) <= TOL
print("  checked: the sweep is non-increasing, and every value is at least")
print(f"  the L-infinity norm, {measures.linf_norm(V)}.")
print()
print("  p = 64 already gives "
      f"{measures.p_norm(V, 64):.9f}, which rounds to 4 at six")
print("  decimal places. The 'limit' arrives quickly in practice.")
print()

print("=" * 74)
print("2. The unit balls: the picture the whole family hangs on")
print("=" * 74)
print()
print("  Every point marked below is at distance 1.0 or less from the centre.")
print("  Same centre, same radius, three different meanings of 'radius'.")
print()
print("    #  inside the L1 ball        -- a DIAMOND")
print("    +  inside L2 but not L1      -- the ring out to the CIRCLE")
print("    .  inside L-inf but not L2   -- the corners of the SQUARE")
print()

WIDTH, HEIGHT = 61, 25
counts = {1: 0, 2: 0, "inf": 0}
for row in range(HEIGHT):
    y = 1.25 - 2.5 * row / (HEIGHT - 1)
    line = []
    for col in range(WIDTH):
        x = -1.25 + 2.5 * col / (WIDTH - 1)
        point = (x, y)
        if measures.p_norm(point, 1) <= 1.0:
            line.append("#")
            counts[1] += 1
            counts[2] += 1
            counts["inf"] += 1
        elif measures.p_norm(point, 2) <= 1.0:
            line.append("+")
            counts[2] += 1
            counts["inf"] += 1
        elif measures.p_norm(point, math.inf) <= 1.0:
            line.append(".")
            counts["inf"] += 1
        else:
            line.append(" ")
    print(("    " + "".join(line)).rstrip())
print()

assert counts[1] < counts[2] < counts["inf"]
print(f"  grid cells inside each ball: L1 {counts[1]}, L2 {counts[2]}, "
      f"L-inf {counts['inf']}")
print("  strictly nested, which is the same fact as the falling column above:")
print("  a bigger p is a more forgiving norm, so its ball is bigger.")
print()
cell_area = (2.5 / (WIDTH - 1)) * (2.5 / (HEIGHT - 1))
print("  Their true areas are 2, pi and 4. Counting cells on this coarse grid")
print(f"  estimates them as {counts[1] * cell_area:.3f}, "
      f"{counts[2] * cell_area:.3f} and {counts['inf'] * cell_area:.3f} --")
print("  close enough to recognise pi, and a good reminder that a picture")
print("  made of characters is an illustration and not a measurement.")
print()
assert abs(counts[2] * cell_area - math.pi) < 0.25

print("=" * 74)
print("3. What has to be true before you may call something a norm")
print("=" * 74)
print()
print("  Four requirements. Any function that satisfies all four is a norm;")
print("  any that misses one is not, whatever it is called.")
print()

v = catalogue.AXIOM_VECTOR
w = catalogue.TRIANGLE_TRIPLE[1]
k = catalogue.AXIOM_SCALAR
zero = (0.0, 0.0, 0.0)
print(f"    v = {v}      w = {w}      k = {k}")
print()

for label, fn in (("L1", measures.l1_norm),
                  ("L2", measures.l2_norm),
                  ("L-infinity", measures.linf_norm)):
    print(f"  {label}")
    nv, nw, nz = fn(v), fn(w), fn(zero)
    print(f"    1. non-negativity        ||v|| = {nv:.6f}  >= 0")
    assert nv >= 0.0 and nw >= 0.0

    print(f"    2. zero only at zero     ||0|| = {nz:.6f}, and ||v|| != 0")
    assert nz == 0.0 and nv > 0.0

    scaled = fn([k * x for x in v])
    print(f"    3. absolute homogeneity  ||{k}v|| = {scaled:.6f}"
          f"  =  |{k}| * ||v|| = {abs(k) * nv:.6f}")
    assert abs(scaled - abs(k) * nv) <= TOL

    summed = fn([a + b for a, b in zip(v, w)])
    print(f"    4. triangle inequality   ||v+w|| = {summed:.6f}"
          f"  <=  ||v|| + ||w|| = {nv + nw:.6f}")
    assert summed <= nv + nw + TOL
    print()

print("  The third is the one people forget. It says doubling a vector must")
print("  exactly double its size -- so 'squared Euclidean distance', which is")
print("  everywhere in machine learning because it avoids a square root, is")
print("  NOT a norm and not a metric. Doubling a vector quadruples it.")
print()
squared = sum(x * x for x in v)
squared_doubled = sum((2 * x) ** 2 for x in v)
print(f"    squared L2 of v      = {squared:.1f}")
print(f"    squared L2 of 2v     = {squared_doubled:.1f}"
      f"   = {squared_doubled / squared:.0f} times, not 2")
assert abs(squared_doubled - 4 * squared) <= TOL
print()
print("  That does not make it useless -- it ranks identically to L2, because")
print("  squaring is monotonic on non-negative numbers, and it is cheaper. It")
print("  makes it useless as a DISTANCE, so never feed it to anything that")
print("  assumes the triangle inequality, such as a ball tree index.")
print()

order_l2 = measures.rank(catalogue.QUERY, catalogue.ARTICLES,
                         measures.l2_distance)
order_sq = measures.rank(
    catalogue.QUERY, catalogue.ARTICLES,
    lambda a, b: sum((x - y) ** 2 for x, y in zip(a, b)))
print(f"    L2 ranking          {[n for n, _ in order_l2]}")
print(f"    squared L2 ranking  {[n for n, _ in order_sq]}")
assert [n for n, _ in order_l2] == [n for n, _ in order_sq]
print("    identical, as promised.")
print()

print("=" * 74)
print("4. numpy.linalg.norm(v, ord=p) IS this family")
print("=" * 74)
print()
print("  `ord` is p. Agreement here is a real check, because measures.py")
print("  computes with `abs`, `**` and `sum` and never calls NumPy.")
print()
print(f"    {'ord':>8}{'measures.py':>16}{'numpy':>16}{'difference':>14}")
print("    " + "-" * 54)
worst = 0.0
for p in (1, 1.5, 2, 3, 8, np.inf):
    mine = measures.p_norm(V, math.inf if np.isinf(p) else p)
    theirs = float(np.linalg.norm(np.asarray(V), ord=p))
    worst = max(worst, abs(mine - theirs))
    label = "inf" if np.isinf(p) else f"{p:g}"
    print(f"    {label:>8}{mine:>16.9f}{theirs:>16.9f}"
          f"{abs(mine - theirs):>14.2e}")
print()
print(f"  worst difference: {worst:.3e}, against a stated tolerance of {TOL:.0e}")
assert worst <= TOL
print()
print("  One difference worth knowing: numpy.linalg.norm accepts ord=0 and")
print("  ord=-1, and neither is a norm. ord=0 counts the non-zero entries,")
print("  which fails absolute homogeneity outright -- doubling a vector does")
print("  not change how many entries are non-zero.")
print()
sparse_values = (0.0, 3.0, 0.0, -7.0)
sparse = np.array(sparse_values)
zero_norm = float(np.linalg.norm(sparse, ord=0))
zero_norm_doubled = float(np.linalg.norm(2 * sparse, ord=0))
print(f"    numpy.linalg.norm({sparse_values}, ord=0) = {zero_norm}")
print(f"    the same vector doubled                        = "
      f"{zero_norm_doubled}")
assert zero_norm == zero_norm_doubled == 2.0
print()
print("  It is still useful -- 'the L0 norm' is how sparsity is counted in")
print("  compressed sensing and in pruning -- and it is still not a norm.")
print("  measures.p_norm refuses p < 1 rather than returning a number:")
print()
try:
    measures.p_norm(V, 0.5)
except ValueError as exc:
    print(f"    p_norm(v, 0.5) -> ValueError: {exc}")
else:  # pragma: no cover - the call above must raise
    raise AssertionError("p_norm should refuse p < 1")
print()

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

"""Which of these measures is a metric, and what it costs you when one is not.

Run from the `examples/` directory:

    ../.venv/bin/python3 03_metrics_and_non_metrics.py

"Metric" is not a compliment. It is a checklist of four properties, and the
reason to care is that indexes, clustering algorithms and proofs of correctness
are built on them. A measure that fails one is not banned -- cosine distance
fails one and is the most used measure in the field -- but you must know which
one it failed and what you were relying on.
"""

from __future__ import annotations

import itertools
import math

import catalogue
import measures
from measures import TOL

print("=" * 74)
print("1. The four axioms")
print("=" * 74)
print()
print("  A function d(x, y) is a METRIC when, for every x, y and z:")
print()
print("    1. non-negativity            d(x, y) >= 0")
print("    2. identity of indiscernibles")
print("                                 d(x, y) = 0 if and only if x = y")
print("    3. symmetry                  d(x, y) = d(y, x)")
print("    4. triangle inequality       d(x, z) <= d(x, y) + d(y, z)")
print()
print("  The fourth is the one with teeth. It says a detour can never be")
print("  shorter than going direct, and it is what lets an index skip whole")
print("  regions of a dataset without looking inside them -- if the query is")
print("  10 away from a cluster centre and the cluster has radius 2, nothing")
print("  in it can be closer than 8, so it need not be opened.")
print()

METRICS = {
    "L1 (Manhattan)": measures.l1_distance,
    "L2 (Euclidean)": measures.l2_distance,
    "L-inf (Chebyshev)": measures.linf_distance,
}

print("=" * 74)
print("2. L1, L2 and L-infinity pass all four")
print("=" * 74)
print()

x, y, z = catalogue.TRIANGLE_TRIPLE
print(f"    x = {x}")
print(f"    y = {y}")
print(f"    z = {z}")
print()

for name, d in METRICS.items():
    print(f"  {name}")
    assert d(x, y) >= 0.0 and d(y, z) >= 0.0 and d(x, z) >= 0.0
    print(f"    1. non-negativity   d(x,y) = {d(x, y):.6f}   >= 0")

    assert d(x, x) == 0.0 and d(x, y) > 0.0
    print(f"    2. zero iff equal   d(x,x) = {d(x, x):.6f}, "
          f"d(x,y) = {d(x, y):.6f}")

    assert abs(d(x, y) - d(y, x)) <= TOL
    print(f"    3. symmetry         d(x,y) = {d(x, y):.6f}"
          f" = d(y,x) = {d(y, x):.6f}")

    # Every one of the six orderings, not just the convenient one.
    worst_slack = math.inf
    for a, b, c in itertools.permutations((x, y, z)):
        slack = d(a, b) + d(b, c) - d(a, c)
        assert slack >= -TOL, (name, slack)
        worst_slack = min(worst_slack, slack)
    print(f"    4. triangle         all 6 orderings hold; tightest slack "
          f"{worst_slack:.6f}")
    print(f"       direct  d(x,z) = {d(x, z):.6f}")
    print(f"       via y            {d(x, y) + d(y, z):.6f}")
    print()

print("=" * 74)
print("3. Cosine distance is NOT a metric, with the counter-example")
print("=" * 74)
print()
print("  Day 103 proved this. Restated here with concrete numbers, because a")
print("  triple you can hold in your head outlasts a proof.")
print()

east, diagonal, north = catalogue.EAST, catalogue.DIAGONAL, catalogue.NORTH
print(f"    east     = {east}      pointing along x")
print(f"    diagonal = {diagonal}      45 degrees between them")
print(f"    north    = {north}      pointing along y")
print()

d_ed = measures.cosine_distance(east, diagonal)
d_dn = measures.cosine_distance(diagonal, north)
d_en = measures.cosine_distance(east, north)
print(f"    cosine_distance(east, diagonal)  = {d_ed:.6f}")
print(f"    cosine_distance(diagonal, north) = {d_dn:.6f}")
print("    ----------------------------------------------")
print(f"    going via the diagonal           = {d_ed + d_dn:.6f}")
print(f"    cosine_distance(east, north)     = {d_en:.6f}   <-- LONGER")
print()
assert d_ed + d_dn < d_en - TOL
print(f"  The direct route is {d_en - (d_ed + d_dn):.6f} longer than the detour.")
print("  No metric may ever allow that. Cosine distance does, so it is a")
print("  DISSIMILARITY and not a distance, whatever the function is called.")
print()
print("  It fails the second axiom too, and this one bites more often:")
print()
doubled = tuple(2 * c for c in east)
print(f"    cosine_distance({east}, {doubled}) = "
      f"{measures.cosine_distance(east, doubled):.6f}")
print(f"    and {east} is not {doubled}")
assert abs(measures.cosine_distance(east, doubled)) <= TOL
assert east != doubled
print()
print("  Distance zero between two things that are not the same thing. For")
print("  cosine that is the FEATURE -- length is what it was asked to ignore --")
print("  but it means cosine cannot tell a document from the same document")
print("  repeated twice, and any deduplication built on it will not either.")
print()
print("  Angular distance, arccos(similarity) / pi, IS a metric on the same")
print("  data and preserves the same ranking, so when an index demands a")
print("  metric that is the standard repair:")
print()
for pair, label in (((east, diagonal), "east / diagonal"),
                    ((diagonal, north), "diagonal / north"),
                    ((east, north), "east / north")):
    ang = math.acos(max(-1.0, min(1.0, measures.cosine_similarity(*pair)))) / math.pi
    print(f"    angular({label:<17}) = {ang:.6f}")
ang_ed = math.acos(measures.cosine_similarity(east, diagonal)) / math.pi
ang_dn = math.acos(measures.cosine_similarity(diagonal, north)) / math.pi
ang_en = math.acos(max(-1.0, min(1.0, measures.cosine_similarity(east, north)))) / math.pi
print(f"    via the diagonal          = {ang_ed + ang_dn:.6f}"
      f"   >=  direct {ang_en:.6f}")
assert ang_ed + ang_dn >= ang_en - TOL
print()

print("=" * 74)
print("4. Jaccard distance and Hamming distance ARE metrics")
print("=" * 74)
print()
print("  Not asserted from a textbook. Checked exhaustively on every triple.")
print()

universe = ("a", "b", "c", "d")
subsets = [frozenset(c)
           for r in range(len(universe) + 1)
           for c in itertools.combinations(universe, r)]
tightest = math.inf
triples = 0
for a, b, c in itertools.product(subsets, repeat=3):
    slack = (measures.jaccard_distance(a, b) + measures.jaccard_distance(b, c)
             - measures.jaccard_distance(a, c))
    assert slack >= -TOL, (a, b, c, slack)
    tightest = min(tightest, slack)
    triples += 1
print(f"  Jaccard distance over all {len(subsets)} subsets of a 4-element set:")
print(f"    {triples} triples checked, none violated the triangle inequality")
print(f"    tightest slack: {tightest:.6f} (0 means equality, which is allowed)")
assert triples == len(subsets) ** 3
print()

strings = list(itertools.product((0, 1), repeat=4))
tightest_h = math.inf
triples_h = 0
for a, b, c in itertools.product(strings, repeat=3):
    slack = (measures.hamming_distance(a, b) + measures.hamming_distance(b, c)
             - measures.hamming_distance(a, c))
    assert slack >= 0, (a, b, c, slack)
    tightest_h = min(tightest_h, slack)
    triples_h += 1
print(f"  Hamming distance over all {len(strings)} 4-bit strings:")
print(f"    {triples_h} triples checked, none violated the triangle inequality")
print(f"    tightest slack: {tightest_h}")
assert triples_h == len(strings) ** 3
print()

print("  The same sweep run on cosine distance finds violations immediately,")
print("  which is what makes the two results above worth having:")
print()
violations = 0
worst = 0.0
vectors = [v for v in itertools.product((0, 1), repeat=4) if any(v)]
for a, b, c in itertools.product(vectors, repeat=3):
    slack = (measures.cosine_distance(a, b) + measures.cosine_distance(b, c)
             - measures.cosine_distance(a, c))
    if slack < -TOL:
        violations += 1
        worst = min(worst, slack)
print(f"    cosine distance over all {len(vectors)} non-zero 4-bit vectors:")
print(f"    {violations} of {len(vectors) ** 3} triples VIOLATE the inequality")
print(f"    worst violation: {worst:.6f}")
assert violations > 0
print()

print("=" * 74)
print("5. What this actually costs")
print("=" * 74)
print()
print("  Metric  ->  ball trees, KD-trees, cover trees, metric-space pruning,")
print("              and the proof that k-medoids terminates.")
print()
print("  Not a metric  ->  none of those are valid. The usual workaround in")
print("              vector databases is to normalise every vector to length 1")
print("              on the way in; once every vector has length 1, cosine")
print("              similarity and Euclidean distance rank identically:")
print()
print("      ||u - v||^2 = 2 - 2 * cosine(u, v)   when ||u|| = ||v|| = 1")
print()
for pair in ((east, diagonal), (diagonal, north), (east, north)):
    u = [c / measures.l2_norm(pair[0]) for c in pair[0]]
    v = [c / measures.l2_norm(pair[1]) for c in pair[1]]
    lhs = measures.l2_distance(u, v) ** 2
    rhs = 2.0 - 2.0 * measures.cosine_similarity(u, v)
    print(f"      {lhs:.9f}  vs  {rhs:.9f}   difference {abs(lhs - rhs):.2e}")
    assert abs(lhs - rhs) <= 1e-9
print()
print("  So the practical advice is not 'avoid cosine'. It is: normalise on")
print("  the way in, then you get cosine's ranking AND a genuine metric, and")
print("  the index is allowed to prune again.")
print()

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

"""Four kinds of data, four right answers. Choose by the shape, not by habit.

Run from the `examples/` directory:

    ../.venv/bin/python3 04_choosing_by_the_shape_of_the_data.py

Grid movement, tolerance checking, categorical fields and sets. Each section
shows a case where the popular default -- Euclidean, or cosine -- is not merely
suboptimal but answers a question nobody asked.
"""

from __future__ import annotations

import catalogue
import measures
from measures import TOL

print("=" * 74)
print("1. Manhattan, Euclidean and Chebyshev on ONE displacement")
print("=" * 74)
print()

a, b = catalogue.FLOOR_FROM, catalogue.FLOOR_TO
dx, dy = b[0] - a[0], b[1] - a[1]
print(f"  A warehouse floor. Go from {a} to {b}, in metres.")
print(f"  The displacement is {dx} across and {dy} along. One pair of points.")
print()
l1 = measures.l1_distance(a, b)
l2 = measures.l2_distance(a, b)
linf = measures.linf_distance(a, b)
print(f"    L1  (Manhattan) = {l1:5.1f}   a picker walking the aisles, one axis")
print("                              at a time. There is no diagonal to walk.")
print(f"    L2  (Euclidean) = {l2:5.1f}   a drone flying it straight. This is the")
print("                              only one that is a physical length here.")
print(f"    Linf (Chebyshev)= {linf:5.1f}   a two-axis gantry whose motors run at")
print("                              the same speed AT THE SAME TIME, so the")
print("                              slower axis alone sets the finishing time.")
print()
assert (l1, l2, linf) == (14.0, 10.0, 8.0)
print("  14, 10 and 8. None is a rounding of another and none is wrong. Each")
print("  is the real cost for a different machine, and if you pick the wrong")
print("  one your route planner optimises a journey nobody takes.")
print()
print(f"  The ordering L-inf <= L2 <= L1 is guaranteed, not a coincidence:")
print("  it is the falling p-column from script 02, applied to a difference.")
assert linf <= l2 <= l1
print()

print("=" * 74)
print("2. Chebyshev, where a single worst component decides alone")
print("=" * 74)
print()
print("  A machined part with four dimensions, in millimetres. It is rejected")
print(f"  if ANY dimension is out by more than {catalogue.PART_TOLERANCE_MM} mm.")
print("  That acceptance rule is an L-infinity ball and cannot be written as")
print("  anything else.")
print()
print(f"    nominal   {catalogue.NOMINAL_PART}")
print()
header = f"    {'batch':<10}{'deviations':<32}{'L1':>8}{'L2':>8}{'L-inf':>8}  verdict"
print(header)
print("    " + "-" * (len(header) - 4))
verdicts = {}
for name, part in catalogue.MEASURED_PARTS.items():
    dev = [round(p - n, 6) for p, n in zip(part, catalogue.NOMINAL_PART)]
    d1 = measures.l1_distance(part, catalogue.NOMINAL_PART)
    d2 = measures.l2_distance(part, catalogue.NOMINAL_PART)
    di = measures.linf_distance(part, catalogue.NOMINAL_PART)
    ok = di <= catalogue.PART_TOLERANCE_MM + TOL
    verdicts[name] = ok
    shown = "[" + ", ".join(f"{v:+.2f}" for v in dev) + "]"
    print(f"    {name:<10}{shown:<32}{d1:>8.2f}{d2:>8.4f}{di:>8.2f}"
          f"  {'ACCEPT' if ok else 'REJECT'}")
print()
assert verdicts == {"batch-A": True, "batch-B": False}
print("  Read that twice. batch-A is out on all four dimensions and its total")
print("  error is nearly double batch-B's -- and batch-A is the one that")
print("  passes. Both L1 and L2 rank batch-B as the better part. Both are")
print("  answering a question the inspection department did not ask.")
print()
print("  Whenever the rule is 'no single feature may be worse than X',")
print("  the measure is Chebyshev. Averaging is not a safe default there;")
print("  it is a way of hiding one bad value behind three good ones.")
print()

print("=" * 74)
print("3. Hamming, for data with no arithmetic in it")
print("=" * 74)
print()
print("  Six categorical fields from a parts register. There is no sense in")
print("  which brass is nearer to steel than nylon is, and any measure that")
print("  subtracts one from the other has invented information.")
print()
print((f"    {'field':<10}"
       + "".join(f"{f:<12}" for f in catalogue.FIELDS)).rstrip())
print((f"    {'reference':<10}"
       + "".join(f"{v:<12}" for v in catalogue.REFERENCE_RECORD)).rstrip())
print()
for name, record in catalogue.CANDIDATE_RECORDS.items():
    marks = "".join(
        f"{(v + ' *') if v != r else v:<12}"
        for v, r in zip(record, catalogue.REFERENCE_RECORD))
    print(f"    {name:<10}{marks}".rstrip())
print()
print("    (* marks a field that differs)")
print()
print(f"    {'record':<10}{'Hamming':>10}{'normalised':>14}")
print("    " + "-" * 30)
hammings = {}
for name, record in catalogue.CANDIDATE_RECORDS.items():
    h = measures.hamming_distance(catalogue.REFERENCE_RECORD, record)
    hammings[name] = h
    print(f"    {name:<10}{h:>10}"
          f"{measures.normalised_hamming(catalogue.REFERENCE_RECORD, record):>14.4f}")
print()
assert hammings == {"part-71": 1, "part-72": 3, "part-73": 6}
print("  part-71 differs only in colour: order it. part-73 shares nothing")
print("  with the reference at all, and the number 6 says exactly that.")
print()
print("  On bits, which is where Hamming defined it in 1950 for error-")
print("  detecting codes, the same count is the number of bit flips between")
print("  two words:")
print()
print(f"    A = {''.join(str(b) for b in catalogue.FLAGS_A)}")
print(f"    B = {''.join(str(b) for b in catalogue.FLAGS_B)}")
diff_marks = "".join(
    "^" if x != y else " "
    for x, y in zip(catalogue.FLAGS_A, catalogue.FLAGS_B))
print(f"        {diff_marks}".rstrip())
flag_h = measures.hamming_distance(catalogue.FLAGS_A, catalogue.FLAGS_B)
print(f"    Hamming distance = {flag_h}")
assert flag_h == 2
print()
print("  And on bits only, Hamming coincides exactly with squared Euclidean")
print("  and with L1, because every difference is 0 or 1 and 1 squared is 1:")
print()
l1_flags = measures.l1_distance(catalogue.FLAGS_A, catalogue.FLAGS_B)
sq_flags = measures.l2_distance(catalogue.FLAGS_A, catalogue.FLAGS_B) ** 2
print(f"    L1 = {l1_flags:.1f}   squared L2 = {sq_flags:.1f}   Hamming = {flag_h}")
assert abs(l1_flags - flag_h) <= TOL and abs(sq_flags - flag_h) <= 1e-9
print()
print("  That coincidence is worth knowing and worth distrusting. It holds")
print("  for BINARY features and collapses the moment a categorical field is")
print("  encoded as 0, 1, 2 -- because then 'nylon' minus 'steel' becomes 2")
print("  and 'brass' minus 'steel' becomes 1, and you have quietly asserted")
print("  that brass is twice as similar to steel as nylon is.")
print()
mislabelled = {"steel": 0, "brass": 1, "nylon": 2}
ref_code = mislabelled[catalogue.REFERENCE_RECORD[0]]
for name in ("part-72", "part-73"):
    code = mislabelled[catalogue.CANDIDATE_RECORDS[name][0]]
    print(f"    integer-encoded material distance, reference to {name}: "
          f"{abs(code - ref_code)}")
print("    ... which is a claim about metallurgy that nobody made.")
print()

print("=" * 74)
print("4. Jaccard against cosine on the same set data")
print("=" * 74)
print()
print("  This is the one most people get wrong, because cosine is the habit.")
print()
q = catalogue.RECIPE_QUERY
print(f"  You want a recipe using: {sorted(q)}")
print()
axes = measures.vocabulary(q, *catalogue.RECIPES.values())
qv = measures.to_binary_vector(q, axes)
print(f"    {'recipe':<14}{'size':>6}{'shared':>8}{'union':>7}"
      f"{'Jaccard':>10}{'cosine':>10}")
print("    " + "-" * 55)
jac, cos = {}, {}
for name, items in catalogue.RECIPES.items():
    jac[name] = measures.jaccard_similarity(q, items)
    cos[name] = measures.cosine_similarity(
        qv, measures.to_binary_vector(items, axes))
    print(f"    {name:<14}{len(items):>6}{len(q & items):>8}"
          f"{len(q | items):>7}{jac[name]:>10.4f}{cos[name]:>10.4f}")
print()

jac_winner = max(jac, key=jac.get)
cos_winner = max(cos, key=cos.get)
print(f"    Jaccard picks {jac_winner}")
print(f"    cosine  picks {cos_winner}")
assert jac_winner == "Shortbread"
assert cos_winner == "Sachertorte"
assert jac_winner != cos_winner
print()
print("  Same two sets. Same query. Opposite answers, and both defensible.")
print()
print("    cosine  = shared / sqrt(|query| * |recipe|)"
      f"  = 4 / sqrt(4*11) = {cos['Sachertorte']:.4f}")
print("    Jaccard = shared / |union|"
      f"                   = 4 / 11        = {jac['Sachertorte']:.4f}")
print()
print("  Sachertorte contains every ingredient you named. Cosine rewards that")
print("  and charges only a square root for the seven extras. Jaccard puts")
print("  the extras in the denominator at full price, so an eleven-ingredient")
print("  cake is not a close match to a four-ingredient request even when it")
print("  is a superset of it.")
print()
print("  Which is right depends on the question:")
print("    'has it got what I asked for?'      -> cosine")
print("    'is it about the same size job?'    -> Jaccard")
print()
print("  For duplicate detection, overlapping tag sets, shingled documents and")
print("  anything where a long item must not out-rank a focused one, Jaccard")
print("  is the safer default -- and unlike cosine distance, 1 - Jaccard is a")
print("  genuine metric, which script 03 checked on all 4096 triples.")
print()
print("  One more asymmetry worth seeing. Cosine on binary data cannot fall")
print("  below Jaccard, ever, because sqrt(|a| * |b|) <= |a union b|:")
print()
pairs_checked = 0
for name, items in catalogue.RECIPES.items():
    assert cos[name] >= jac[name] - TOL
    pairs_checked += 1
    print(f"    {name:<14} cosine {cos[name]:.4f}  >=  Jaccard {jac[name]:.4f}")
assert pairs_checked == 2
print()
print("  So cosine is systematically the more generous of the two on sets.")
print("  If your relevance scores look suspiciously high, that is a candidate")
print("  explanation before you go looking for a bug.")
print()

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

"""Mahalanobis: Euclidean distance after accounting for how the data varies.

Run from the `examples/` directory:

    ../.venv/bin/python3 05_mahalanobis_distance.py

Euclidean distance treats every direction as equally surprising. Real data
does not: it has a grain, and moving along the grain is ordinary while moving
across it is an event. Mahalanobis distance is what you get when you measure in
the data's own units instead of the axes' units, and Day 106's eigenvectors of
the covariance matrix are exactly the directions it measures along.
"""

from __future__ import annotations

import math

import numpy as np

import catalogue
import measures
from measures import TOL

DATA = catalogue.SENSOR_READINGS

print("=" * 74)
print("1. Eight readings from two sensors that move together")
print("=" * 74)
print()
print("    reading   sensor A   sensor B")
print("    " + "-" * 30)
for i, (sa, sb) in enumerate(DATA, start=1):
    print(f"    {i:<10}{sa:>9.1f}{sb:>11.1f}")
print()

mean = measures.column_means(DATA)
print(f"  mean = {tuple(mean)}")
assert mean == [0.0, 0.0]
print()

# A picture of the grain, drawn from the data rather than described.
print("  Plotted, with the two probe points marked:")
print()
LO, HI = -5, 5
for row in range(HI, LO - 1, -1):
    line = []
    for col in range(LO, HI + 1):
        point = (float(col), float(row))
        if point == catalogue.PROBE_ALONG:
            line.append(" A")
        elif point == catalogue.PROBE_ACROSS:
            line.append(" X")
        elif point in DATA:
            line.append(" o")
        elif col == 0 and row == 0:
            line.append(" +")
        elif col == 0:
            line.append(" |")
        elif row == 0:
            line.append(" -")
        else:
            line.append(" .")
    print("      " + "".join(line))
print()
print("      o  a reading         +  the mean")
print(f"      A  probe {catalogue.PROBE_ALONG}, ALONG the grain")
print(f"      X  probe {catalogue.PROBE_ACROSS}, ACROSS it")
print()
print("  Every reading sits near the line B = A. The two sensors agree, all")
print("  day, and the eight points say so without anyone writing it down.")
print()

print("=" * 74)
print("2. The covariance matrix")
print("=" * 74)
print()
cov = measures.covariance_matrix(DATA)
print("    covariance = [[%.4f, %.4f]," % (cov[0][0], cov[0][1]))
print("                  [%.4f, %.4f]]" % (cov[1][0], cov[1][1]))
print()
assert cov == [[7.5, 7.0], [7.0, 7.5]]
print("  Exactly [[7.5, 7.0], [7.0, 7.5]], with no floating-point residue,")
print("  because the data was chosen so you can check it by hand:")
print()
print("    variance of A  = (16+9+4+1+1+4+9+16) / 8 = 60 / 8 = 7.5")
print("    covariance A,B = (12+12+2+2+2+2+12+12) / 8 = 56 / 8 = 7.0")
print()
corr = cov[0][1] / math.sqrt(cov[0][0] * cov[1][1])
print(f"  correlation = 7.0 / 7.5 = {corr:.6f}   -- very nearly 1")
assert abs(corr - 7.0 / 7.5) <= TOL
print()

cov_np = np.cov(np.asarray(DATA, dtype=float), rowvar=False, bias=True)
print("  NumPy computes the same matrix. `bias=True` is the population")
print("  divisor n, which is what measures.covariance_matrix uses and what")
print("  scikit-learn's StandardScaler uses; the default `bias=False` divides")
print("  by n - 1 and is a different, also-correct, answer to a different")
print("  question.")
print()
print(f"    numpy (bias=True)  = {cov_np.tolist()}")
print(f"    numpy (bias=False) = "
      f"{np.cov(np.asarray(DATA, dtype=float), rowvar=False).tolist()}")
assert np.allclose(cov_np, np.asarray(cov), atol=TOL)
print()

print("=" * 74)
print("3. The inverse, in pure Python, checked against NumPy")
print("=" * 74)
print()
inv = measures.inverse(cov)
det = cov[0][0] * cov[1][1] - cov[0][1] * cov[1][0]
print(f"    determinant = 7.5*7.5 - 7.0*7.0 = {det}")
print("    inverse     = [[%.6f, %.6f]," % (inv[0][0], inv[0][1]))
print("                   [%.6f, %.6f]]" % (inv[1][0], inv[1][1]))
print()
inv_np = np.linalg.inv(np.asarray(cov))
worst = float(np.max(np.abs(np.asarray(inv) - inv_np)))
print(f"    numpy.linalg.inv agrees to {worst:.3e}, tolerance {TOL:.0e}")
assert worst <= TOL
print()
print("  measures.inverse is Gauss-Jordan elimination written out by hand, so")
print("  the Mahalanobis numbers below owe nothing to NumPy and agreeing with")
print("  NumPy means something.")
print()
product = measures.matmul(cov, inv)
print(f"    covariance * inverse = {[[round(x, 12) for x in r] for r in product]}")
assert abs(product[0][0] - 1) <= TOL and abs(product[1][1] - 1) <= TOL
assert abs(product[0][1]) <= TOL and abs(product[1][0]) <= TOL
print()

print("=" * 74)
print("4. Two points Euclidean cannot tell apart")
print("=" * 74)
print()
along, across = catalogue.PROBE_ALONG, catalogue.PROBE_ACROSS
eu_along = measures.l2_distance(along, mean)
eu_across = measures.l2_distance(across, mean)
ma_along = measures.mahalanobis_distance(along, mean, inv)
ma_across = measures.mahalanobis_distance(across, mean, inv)

print(f"    {'probe':<14}{'Euclidean':>12}{'Mahalanobis':>14}")
print("    " + "-" * 40)
print(f"    {str(along):<14}{eu_along:>12.6f}{ma_along:>14.6f}")
print(f"    {str(across):<14}{eu_across:>12.6f}{ma_across:>14.6f}")
print()
assert abs(eu_along - eu_across) <= TOL
assert abs(eu_along - math.sqrt(18.0)) <= TOL
print(f"  Euclidean: identical, both sqrt(18) = {eu_along:.6f}. It has no way")
print("  to distinguish them, because it does not know the data exists.")
print()
assert ma_across > ma_along
print(f"  Mahalanobis: {ma_along:.6f} against {ma_across:.6f}, a factor of "
      f"{ma_across / ma_along:.4f}.")
print("  Both sensors reading 3 together is a perfectly ordinary Tuesday. One")
print("  reading +3 while the other reads -3 has never happened in this")
print("  dataset, and the number says so.")
print()
print("  The second value is 6, and here is why every comparison in this lab")
print("  states a tolerance. Two correct implementations of the same inverse")
print("  disagree in the last bit, and it survives all the way to the answer:")
print()
ma_np = float(math.sqrt(
    np.asarray(across) @ inv_np @ np.asarray(across)))
print(f"    via measures.inverse   (Gauss-Jordan)  {ma_across!r}")
print(f"    via numpy.linalg.inv   (LAPACK)        {ma_np!r}")
print(f"    difference                             {abs(ma_across - ma_np):.3e}")
assert abs(ma_across - 6.0) <= TOL
assert abs(ma_np - 6.0) <= TOL
print()
print("  Neither is wrong and neither is 'more accurate'. `== 6.0` would pass")
print("  for one and fail for the other, which is the entire argument against")
print("  writing `==` between two floats you did not personally construct.")
print()
print("  That is the whole argument for the measure. An anomaly detector built")
print("  on Euclidean distance has to score these two the same. One of them is")
print("  a sensor fault.")
print()

print("=" * 74)
print("5. Where the numbers come from: Day 106's eigenvectors")
print("=" * 74)
print()
values, vectors = np.linalg.eigh(np.asarray(cov))
print("  Eigen-decomposition of the covariance matrix:")
for value, vector in zip(values, vectors.T):
    print(f"    eigenvalue {value:8.4f}   eigenvector "
          f"({vector[0]:+.6f}, {vector[1]:+.6f})")
print()
assert abs(sorted(values)[0] - 0.5) <= 1e-12
assert abs(sorted(values)[1] - 14.5) <= 1e-12
print("  0.5 and 14.5. The large one belongs to the (1, 1) direction -- along")
print("  the grain, where the data spreads a lot -- and the small one to")
print("  (1, -1), across it, where the data barely spreads at all. An")
print("  eigenvector's SIGN is arbitrary, which is why NumPy prints the small")
print("  one as (-0.707107, +0.707107): that is the same line as (1, -1),")
print("  pointing the other way, and no distance below can tell the")
print("  difference because every component is squared.")
print()
print("  Mahalanobis distance is Euclidean distance measured in those")
print("  directions, with each one divided by the square root of its own")
print("  eigenvalue. Worked by hand for both probes:")
print()
axis_along = (1 / math.sqrt(2), 1 / math.sqrt(2))
axis_across = (1 / math.sqrt(2), -1 / math.sqrt(2))
for label, probe in (("along  (3, 3)", along), ("across (3, -3)", across)):
    c_big = measures.dot(probe, axis_along)
    c_small = measures.dot(probe, axis_across)
    by_hand = math.sqrt(c_big ** 2 / 14.5 + c_small ** 2 / 0.5)
    direct = measures.mahalanobis_distance(probe, mean, inv)
    print(f"    {label}")
    print(f"      component along  (1, 1)/sqrt(2) = {c_big:+.6f}"
          f"   / sqrt(14.5) = {c_big / math.sqrt(14.5):+.6f}")
    print(f"      component across (1,-1)/sqrt(2) = {c_small:+.6f}"
          f"   / sqrt( 0.5) = {c_small / math.sqrt(0.5):+.6f}")
    print(f"      hypotenuse of those two         = {by_hand:.6f}")
    print(f"      mahalanobis_distance says       = {direct:.6f}")
    assert abs(by_hand - direct) <= 1e-9
    print()
print("  So Mahalanobis is not a new kind of distance at all. It is Euclidean")
print("  distance in a coordinate system the data chose for itself, and the")
print("  eigenvectors Day 106 built are the axes of that system.")
print()

print("=" * 74)
print("6. Substituting the identity gives back Euclidean, exactly")
print("=" * 74)
print()
identity = [[1.0, 0.0], [0.0, 1.0]]
worst = 0.0
for probe in (along, across, (1.0, 0.0), (-2.5, 4.75), (0.0, 0.0)):
    a_val = measures.mahalanobis_distance(probe, mean, identity)
    b_val = measures.l2_distance(probe, mean)
    worst = max(worst, abs(a_val - b_val))
    print(f"    probe {str(probe):<14} mahalanobis {a_val:>10.6f}"
          f"   euclidean {b_val:>10.6f}")
print()
print(f"  worst difference: {worst:.3e}")
assert worst <= TOL
print()
print("  Which is the cleanest way to see what the covariance is doing: it is")
print("  the thing that would be the identity if every feature had variance 1")
print("  and no feature had anything to do with any other. Real data is never")
print("  that, and Euclidean distance quietly assumes it always is.")
print()
print("  The cost is real and worth stating. Mahalanobis needs an invertible")
print("  covariance matrix, which needs more rows than columns and no two")
print("  features that are exact duplicates -- and it needs re-estimating when")
print("  the data drifts. A singular covariance raises here rather than")
print("  returning a plausible number:")
print()
duplicate_feature = [(1.0, 2.0), (2.0, 4.0), (3.0, 6.0), (4.0, 8.0)]
try:
    measures.inverse(measures.covariance_matrix(duplicate_feature))
except ValueError as exc:
    print(f"    second feature = 2 * first  ->  ValueError: {exc}")
else:  # pragma: no cover - the call above must raise
    raise AssertionError("a singular covariance should refuse to invert")
print()

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

"""The units in your table decide your answer, and nobody voted on them.

Run from the `examples/` directory:

    ../.venv/bin/python3 06_scaling_changes_the_answer.py

Every distance in this lab sums contributions across features. Nothing in that
sum knows that one column is in metres and another in grams, so the column with
the bigger numbers wins the argument -- not by being more important, but by
being written down in smaller units.
"""

from __future__ import annotations

import numpy as np

import catalogue
import measures
from measures import TOL

B = catalogue.BEARINGS
Q = catalogue.BEARING_QUERY
ROWS = list(B.values())

print("=" * 74)
print("1. A bearing catalogue in the units the supplier used")
print("=" * 74)
print()
print(f"    {'part':<8}{catalogue.BEARING_FEATURES[0]:>22}"
      f"{catalogue.BEARING_FEATURES[1]:>14}")
print("    " + "-" * 44)
print(f"    {'WANTED':<8}{Q[0]:>22.3f}{Q[1]:>14.1f}")
print("    " + "-" * 44)
for name, row in B.items():
    print(f"    {name:<8}{row[0]:>22.3f}{row[1]:>14.1f}")
print()
print("  Bore diameter is recorded in METRES, so every number in that column")
print("  is about 0.02. Mass is in GRAMS, so every number in that one is in")
print("  the hundreds. Both columns matter to an engineer. Only one of them")
print("  is going to matter to a Euclidean distance.")
print()

print("=" * 74)
print("2. Rank on the raw numbers")
print("=" * 74)
print()
raw = measures.rank(Q, B, measures.l2_distance)
head2 = (f"    {'part':<8}{'distance':>14}{'bore term':>16}"
         f"{'mass term':>16}{'bore share':>13}")
print(head2)
print("    " + "-" * (len(head2) - 4))
for name, d in raw:
    bore = (Q[0] - B[name][0]) ** 2
    mass = (Q[1] - B[name][1]) ** 2
    share = bore / (bore + mass) if bore + mass else 0.0
    print(f"    {name:<8}{d:>14.6f}{bore:>16.2e}{mass:>16.2f}{share:>12.6%}")
print()
raw_winner = raw[0][0]
print(f"  Winner: {raw_winner}")
assert raw_winner == "R"
print()
print("  Look at the last column. The bore diameter contributes less than one")
print("  ten-thousandth of one per cent of every distance in the table. This")
print("  is not a ranking on two features. It is a ranking on mass, with a")
print("  rounding error attached.")
print()
print(f"  And R is unusable. The query wants a {Q[0] * 1000:.0f} mm bore; R has "
      f"a {B['R'][0] * 1000:.0f} mm bore,")
print("  60 per cent oversize, a part that will not fit the shaft. It wins")
print("  because it is 2 g from the target mass, and mass is the only thing")
print("  being measured.")
print()
print(f"  P has EXACTLY the bore asked for and comes "
      f"{[n for n, _ in raw].index('P') + 1} of {len(raw)}.")
print()

print("=" * 74)
print("3. The same ranking after standardising")
print("=" * 74)
print()
means = measures.column_means(ROWS)
stds = measures.column_stds(ROWS)
print(f"    column means               {[round(m, 6) for m in means]}")
print(f"    column standard deviations {[round(s, 6) for s in stds]}")
print()
print("  The query is standardised with the CATALOGUE's numbers, not its own.")
print("  Standardising a single row against itself gives a row of zeros, which")
print("  is a mistake with a long history in production retrieval systems.")
print()

q_std = measures.standardise([Q], means, stds)[0]
b_std = {name: measures.standardise([row], means, stds)[0]
         for name, row in B.items()}
print(f"    {'part':<8}{'bore (z)':>12}{'mass (z)':>12}{'distance':>14}")
print("    " + "-" * 46)
scaled = measures.rank(q_std, b_std, measures.l2_distance)
print(f"    {'WANTED':<8}{q_std[0]:>12.4f}{q_std[1]:>12.4f}")
for name, d in scaled:
    print(f"    {name:<8}{b_std[name][0]:>12.4f}{b_std[name][1]:>12.4f}"
          f"{d:>14.6f}")
print()
scaled_winner = scaled[0][0]
print(f"  Winner: {scaled_winner}")
assert scaled_winner == "P"
assert scaled_winner != raw_winner
print()
print(f"  The winner changed from {raw_winner} to {scaled_winner}. Same six")
print("  parts, same query, same Euclidean distance, same code. The only")
print("  thing that changed is that both columns now speak in standard")
print("  deviations of the catalogue, so a 12 mm bore error costs what a 12 mm")
print("  bore error is worth rather than what it looks like next to 40 grams.")
print()
print(f"    {'part':<8}{'raw rank':>11}{'standardised rank':>20}{'moved':>8}")
print("    " + "-" * 47)
raw_order = [n for n, _ in raw]
std_order = [n for n, _ in scaled]
moved = 0
for name in B:
    a, c = raw_order.index(name) + 1, std_order.index(name) + 1
    if a != c:
        moved += 1
    print(f"    {name:<8}{a:>11}{c:>20}{('yes' if a != c else '-'):>8}")
print()
print(f"  {moved} of the {len(B)} parts moved -- and they are the two the")
print("  decision is between. P and R swap places, first for third.")
assert moved == 2
assert (raw_order.index("P"), std_order.index("P")) == (2, 0)
assert (raw_order.index("R"), std_order.index("R")) == (0, 2)
print()

print("=" * 74)
print("4. It is the UNITS, not the standardising")
print("=" * 74)
print()
print("  The clearest proof that the raw ranking was an artefact: change no")
print("  data at all, only the unit the bore column is written in.")
print()
for label, factor in (("metres", 1.0), ("millimetres", 1e3),
                      ("micrometres", 1e6)):
    q_u = (Q[0] * factor, Q[1])
    b_u = {n: (v[0] * factor, v[1]) for n, v in B.items()}
    order = [n for n, _ in measures.rank(q_u, b_u, measures.l2_distance)]
    print(f"    bore in {label:<13} {order}")
    if label == "micrometres":
        assert order[0] == "P"
    if label == "metres":
        assert order[0] == "R"
print()
print("  In metres the answer is R. In micrometres the answer is P. The parts")
print("  did not change; a column header did. Any pipeline that does not")
print("  normalise is quietly letting whoever chose the units decide the")
print("  ranking, and that person was usually not thinking about distances.")
print()

print("=" * 74)
print("5. Standardising is not the only choice, and it is not free")
print("=" * 74)
print()
lo = [min(r[j] for r in ROWS) for j in range(len(ROWS[0]))]
hi = [max(r[j] for r in ROWS) for j in range(len(ROWS[0]))]


def min_max(row):
    return [(row[j] - lo[j]) / (hi[j] - lo[j]) for j in range(len(row))]


mm_order = [n for n, _ in measures.rank(
    min_max(Q), {n: min_max(v) for n, v in B.items()}, measures.l2_distance)]
print(f"    z-score (mean 0, sd 1)     {std_order}")
print(f"    min-max (squashed to 0-1) {mm_order}")
assert mm_order[0] == "P"
print()
print("  Both agree here, which will not always happen. The trade-off:")
print()
print("    z-score   assumes nothing about the range, so an outlier stretches")
print("              the standard deviation and squashes everything else")
print("              toward zero. Handles unbounded features.")
print()
print("    min-max   pins the range to 0-1 exactly, which is what an image")
print("              pipeline usually wants -- and one outlier now decides the")
print("              WHOLE scale, and a value outside the training range comes")
print("              out above 1 or below 0.")
print()
print("  There is a third answer that people forget: do not scale, and choose")
print("  a measure that does not need it. Mahalanobis divides by the data's")
print("  own spread as part of its definition, so it needs no scaling step at")
print("  all -- and it does NOT give the same answer, which is worth more than")
print("  if it had:")
print()
cov_inv = measures.inverse(measures.covariance_matrix(ROWS))
maha_order = [n for n, _ in measures.rank(
    Q, B, lambda a, b: measures.mahalanobis_distance(a, b, cov_inv))]
print(f"    Mahalanobis on the RAW numbers  {maha_order}")
print(f"    z-score then Euclidean          {std_order}")
print(f"    raw Euclidean                   {raw_order}")
assert maha_order[0] == "U"
assert maha_order[1] == "P"
assert maha_order.index("R") == 4
print()
print("  Both cures demote the unusable part: R falls from 1st to 3rd under")
print("  standardising and to 5th under Mahalanobis. They disagree at the top,")
print("  where Mahalanobis prefers U and standardising prefers P.")
print()
print("  The disagreement is not noise, and it is the reason to know both.")
print("  In this catalogue bore and mass are correlated -- bigger bearings are")
print("  heavier -- and Mahalanobis removes that shared movement before")
print("  measuring, while standardising only rescales each column separately.")
print()
cov_raw = measures.covariance_matrix(ROWS)
corr = cov_raw[0][1] / (stds[0] * stds[1])
print(f"    correlation between bore and mass: {corr:+.4f}")
assert corr > 0.7
print()
print("  P is close to the query in bore and 40 g heavy. Once you know that")
print("  heavier goes with wider in this catalogue, being wide-for-its-mass or")
print("  heavy-for-its-bore is the surprising thing, and U -- which is smaller")
print("  and lighter TOGETHER, along the grain -- reads as the nearer part.")
print("  Whether you want that is a modelling decision, which is the day's")
print("  entire subject.")
print()

print("=" * 74)
print("6. Is this catalogue cherry-picked? A seeded sweep says no")
print("=" * 74)
print()
print("  Six parts and one query were chosen by hand to make the point")
print("  legible. Here is the same experiment on random catalogues, drawn")
print("  from numpy.random.default_rng(107) -- a SEEDED generator, so this")
print("  run reproduces on this machine, and the claim asserted below is a")
print("  RANGE rather than an exact count, because NumPy does not promise")
print("  that a generator's stream survives a version change.")
print()
rng = np.random.default_rng(107)
TRIALS = 2000
spread = np.array([0.04, 500.0])
flips = 0
for _ in range(TRIALS):
    cat = rng.random((6, 2)) * spread
    query = rng.random(2) * spread
    raw_best = int(np.argmin(np.linalg.norm(cat - query, axis=1)))
    mu, sd = cat.mean(axis=0), cat.std(axis=0)
    std_best = int(np.argmin(
        np.linalg.norm((cat - mu) / sd - (query - mu) / sd, axis=1)))
    flips += raw_best != std_best
print(f"    {TRIALS} random catalogues, same two units")
print(f"    the winner changed after standardising in {flips} of them"
      f"  ({flips / TRIALS:.1%})")
print()
assert 0.35 <= flips / TRIALS <= 0.75, flips
print("  Between a third and three quarters, every time this has been run.")
print("  Standardising is not a tweak that occasionally matters. On features")
print("  in mismatched units it decides the answer about half the time.")
print()

print("=" * 74)
print("7. Where cosine sits in this")
print("=" * 74)
print()
print("  Cosine similarity is often described as 'scale invariant', and that")
print("  is true of the wrong scale. It ignores the length of a VECTOR. It")
print("  does not ignore the units of a COLUMN, and it cannot, because")
print("  changing one column's units rotates every vector in the table.")
print()
for label, factor in (("metres", 1.0), ("micrometres", 1e6)):
    q_u = (Q[0] * factor, Q[1])
    b_u = {n: (v[0] * factor, v[1]) for n, v in B.items()}
    order = [n for n, _ in measures.rank(q_u, b_u, measures.cosine_similarity,
                                         higher_is_better=True)]
    print(f"    cosine, bore in {label:<13} {order}")
    if label == "metres":
        cosine_metres = order
    else:
        cosine_micro = order
assert cosine_metres != cosine_micro
print()
print("  Different order, same data. 'Scale invariant' is a claim about")
print("  multiplying a whole vector by a constant, and it is worth knowing")
print("  exactly that much and no more.")
print()
doubled = {n: tuple(2 * c for c in v) for n, v in B.items()}
same = measures.rank(Q, doubled, measures.cosine_similarity,
                     higher_is_better=True)
assert [n for n, _ in same] == cosine_metres
print("  What it IS invariant to, checked: doubling every candidate vector")
print("  leaves the cosine ranking untouched.")
print()
worst = max(abs(measures.cosine_similarity(Q, doubled[n])
                - measures.cosine_similarity(Q, B[n])) for n in B)
print(f"    largest change in any cosine score: {worst:.3e}"
      f"   (tolerance {TOL:.0e})")
assert worst <= TOL
print()

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

"""The data this lab argues over. Written by hand, on purpose.

Nothing here is downloaded and nothing here is random. Every number is small
enough to check on paper, and each dataset was chosen because it makes exactly
one measure look right and the others look wrong -- which is the point of the
day.

There IS a seeded random generator in this lab, in
`06_scaling_changes_the_answer.py`, used to show that the scaling effect is not
a property of these six hand-picked parts. It is seeded with
`numpy.random.default_rng(107)` and every claim made about it is structural (a
count, a direction, a percentage floor) rather than a specific digit, because
NumPy does not promise that a generator's exact stream survives a version
change. Everything asserted to the last decimal place in this lab comes from the
literal tables below.
"""

from __future__ import annotations

# ---------------------------------------------------------------------------
# 1. The opening disagreement: three articles, one query, three winners.
# ---------------------------------------------------------------------------
#
# A tiny help-centre search. Each article is counted over four terms, and the
# query is the reader's own short note. Raw counts, not frequencies -- which is
# exactly what makes the three measures disagree.

TERMS = ("norm", "distance", "vector", "cluster")

QUERY = (4, 3, 2, 1)

# "Aisle" mentions cluster six times where the query mentions it once, and
# matches the other three terms exactly. One big disagreement, nothing else.
#
# "Beacon" is a little off on three of the four terms and exact on the fourth.
# Three small disagreements adding to more in total than Aisle's single one.
#
# "Cartogram" is the query's profile at exactly three times the length: a long
# article on precisely this topic. Its direction is identical, so its cosine
# similarity is exactly 1.0, and its raw counts are further away than anything
# else here.
ARTICLES: dict[str, tuple[int, ...]] = {
    "Aisle": (4, 3, 2, 6),
    "Beacon": (6, 1, 4, 1),
    "Cartogram": (12, 9, 6, 3),
}

# ---------------------------------------------------------------------------
# 2. Chebyshev, and where a single worst component decides.
# ---------------------------------------------------------------------------
#
# One displacement, in metres, across a warehouse floor laid out on aisles.
# The same two points, three operationally different answers:
#
#   L1  = 14  a picker who must walk the aisles, one axis at a time
#   L2  = 10  a drone that can fly the diagonal
#   Linf = 8  a two-axis gantry whose motors run at once, so the slower axis
#             alone sets the finishing time
#
# None of the three is the "real" distance. Each is the real distance for a
# different machine.

FLOOR_FROM = (0.0, 0.0)
FLOOR_TO = (6.0, 8.0)

# A machined part and its nominal dimensions, in millimetres. The part is
# rejected if ANY dimension is out by more than the tolerance -- which is a
# Chebyshev ball, and nothing else.
NOMINAL_PART = (40.00, 25.00, 12.00, 6.00)
PART_TOLERANCE_MM = 0.05
MEASURED_PARTS: dict[str, tuple[float, ...]] = {
    # Four dimensions each a little out. Total error is large, worst is small.
    "batch-A": (40.04, 24.96, 12.04, 5.96),
    # Three dimensions perfect, one badly out. Total error is smaller.
    "batch-B": (40.00, 25.00, 12.00, 6.09),
}

# ---------------------------------------------------------------------------
# 3. Hamming, for data with no arithmetic in it.
# ---------------------------------------------------------------------------
#
# Six categorical fields from a parts register. Subtracting "steel" from
# "brass" is not a smaller number than subtracting "steel" from "nylon"; it is
# not a number at all. Hamming counts the fields that differ and refuses to
# invent an ordering.

FIELDS = ("material", "finish", "thread", "grade", "colour", "origin")

REFERENCE_RECORD = ("steel", "zinc", "M8", "8.8", "silver", "IN")
CANDIDATE_RECORDS: dict[str, tuple[str, ...]] = {
    "part-71": ("steel", "zinc", "M8", "8.8", "black", "IN"),
    "part-72": ("brass", "zinc", "M8", "10.9", "silver", "DE"),
    "part-73": ("nylon", "plain", "M6", "4.6", "white", "CN"),
}

# The same measure on bits, which is where Hamming was defined: two 8-bit
# feature flags from the same register.
FLAGS_A = (1, 0, 1, 1, 0, 0, 1, 0)
FLAGS_B = (1, 0, 0, 1, 0, 1, 1, 0)

# ---------------------------------------------------------------------------
# 4. Jaccard against cosine on the same set-like data.
# ---------------------------------------------------------------------------
#
# Ingredient lists. The query has four ingredients.
#
#   "Sachertorte" contains ALL FOUR of them, plus seven more.
#   "Shortbread"  shares two of the four and has one ingredient of its own.
#
# Cosine, which divides by the square root of the two sizes, prefers
# Sachertorte: everything asked for is present. Jaccard, which counts the union
# in the denominator, prefers Shortbread: Sachertorte's seven extra ingredients
# are seven things the two recipes do not share, and Jaccard charges for them.
#
# Neither is wrong. They answer different questions -- "is what I asked for
# there?" against "how much of everything involved is shared?" -- and the day's
# job is to notice that you have to pick one.

RECIPE_QUERY = frozenset({"flour", "butter", "sugar", "egg"})
RECIPES: dict[str, frozenset[str]] = {
    "Sachertorte": frozenset({
        "flour", "butter", "sugar", "egg",
        "cocoa", "apricot jam", "chocolate", "vanilla",
        "salt", "milk", "almond",
    }),
    "Shortbread": frozenset({"flour", "butter", "cornflour"}),
}

# ---------------------------------------------------------------------------
# 5. Mahalanobis: Euclidean after accounting for how the data actually varies.
# ---------------------------------------------------------------------------
#
# Eight readings of two sensors that move together almost perfectly. The mean
# is exactly (0, 0) and the population covariance comes out exactly
#
#     [[7.5, 7.0],
#      [7.0, 7.5]]
#
# whose determinant is exactly 7.25. The data lies along the line y = x: that
# is the grain of it, and Day 106's eigenvectors of this matrix are what name
# that direction.
#
# The two probe points are the same Euclidean distance from the mean -- both
# sqrt(18) = 4.2426... -- and nothing about ordinary distance can tell them
# apart. Mahalanobis can: ALONG the grain is cheap, ACROSS it is expensive.

SENSOR_READINGS: tuple[tuple[float, float], ...] = (
    (-4.0, -3.0),
    (-3.0, -4.0),
    (-2.0, -1.0),
    (-1.0, -2.0),
    (1.0, 2.0),
    (2.0, 1.0),
    (3.0, 4.0),
    (4.0, 3.0),
)

# Along the grain of the data: both sensors high together, which is what this
# pair of sensors does all day.
PROBE_ALONG = (3.0, 3.0)
# Across the grain: one sensor high while the other is low, which never happens
# in the eight readings above. Same Euclidean distance. Not the same event.
PROBE_ACROSS = (3.0, -3.0)

# ---------------------------------------------------------------------------
# 6. The scaling demonstration: the thing that silently decides your answer.
# ---------------------------------------------------------------------------
#
# A bearing catalogue with two features recorded in the units the supplier
# happened to use: bore diameter in METRES and mass in GRAMS. The numbers in
# one column are around 0.02 and in the other around 350, so squared
# differences in the second column are roughly ten million times larger. The
# bore column does not lose the argument. It never enters it.

BEARING_FEATURES = ("bore diameter (m)", "mass (g)")

BEARING_QUERY = (0.020, 300.0)

BEARINGS: dict[str, tuple[float, float]] = {
    # Bore matches the query EXACTLY. 40 g heavier.
    "P": (0.020, 340.0),
    # Bore is 12 mm too big -- 60 per cent out, and unusable. Mass is 2 g off.
    "R": (0.032, 302.0),
    "S": (0.008, 250.0),
    "T": (0.026, 410.0),
    "U": (0.014, 275.0),
    "V": (0.038, 500.0),
}

# ---------------------------------------------------------------------------
# 7. The counter-example that shows cosine distance is not a metric.
# ---------------------------------------------------------------------------
#
# Day 103 proved this. It is restated here rather than re-derived, because a
# concrete triple is worth more than the proof once you have seen the proof.
#
#   cosine_distance(EAST, DIAGONAL) + cosine_distance(DIAGONAL, NORTH)
#     = 0.2929 + 0.2929 = 0.5858
#   cosine_distance(EAST, NORTH)
#     = 1.0
#
# The direct route is longer than going via a third point, which no metric may
# ever allow.

EAST = (1.0, 0.0)
DIAGONAL = (1.0, 1.0)
NORTH = (0.0, 1.0)

# The triple used for the POSITIVE side of the same check: L1, L2 and
# L-infinity must all satisfy the triangle inequality on every triple, and
# these three vectors are checked exhaustively in all six orderings.
TRIANGLE_TRIPLE = ((1.0, 7.0, 2.0), (4.0, 1.0, 9.0), (-2.0, 3.0, 3.0))

# The single vector every norm axiom is checked on, and the scalar it is
# multiplied by for absolute homogeneity.
AXIOM_VECTOR = (3.0, -4.0, 12.0)
AXIOM_SCALAR = -2.5

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

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

"""The reference implementation: every measure in this lab, in pure Python.

Read this AFTER you have attempted `starter/measures.py`. Nothing here uses
NumPy. That is deliberate and it is the whole basis of the day's evidence: if
these functions were built out of NumPy calls, then checking them against
NumPy would be checking NumPy against itself and would prove nothing at all.

NumPy appears in the tests and in the demonstration scripts, where it is the
independent answer. `numpy.linalg.norm(v, ord=p)` is exactly the p-norm family
implemented below, and agreeing with it to 1e-12 on values this code computed
from `abs`, `**` and `sum` is a real check.

Two conventions, fixed and used everywhere:

1. A vector is a plain sequence of floats. Two vectors compared must have the
   same length, and every function here says so by raising rather than by
   zipping the shorter one and quietly answering the wrong question.

2. A DISTANCE gets smaller as things get more alike; a SIMILARITY gets larger.
   Every function name says which it is. Mixing them up is the single most
   common way to build a retrieval system that returns the worst match with
   great confidence, so this module never guesses: `rank` takes an explicit
   `higher_is_better` flag.
"""

from __future__ import annotations

import math
from collections.abc import Callable, Hashable, Iterable, Sequence

Vector = Sequence[float]
Matrix = list[list[float]]

# Every float comparison in this lab is made against this tolerance. It is not
# decoration. `math.sqrt(2) ** 2` is 2.0000000000000004, and the Mahalanobis
# result this lab is built around comes out as exactly 6.0 through the
# Gauss-Jordan inverse below and as 5.999999999999999 through
# `numpy.linalg.inv` -- two correct routes to the same number, disagreeing in
# the last bit. `== 6.0` would pass for one and fail for the other.
TOL = 1e-12


class DimensionMismatch(ValueError):
    """Raised when two vectors of different lengths are compared.

    Subclasses ValueError so that an existing `except ValueError` catches it,
    matching how NumPy reports the same mistake.
    """


def _paired(u: Vector, v: Vector) -> list[tuple[float, float]]:
    """Pair two vectors elementwise, refusing to compare different lengths."""
    if len(u) != len(v):
        raise DimensionMismatch(
            f"vectors have different lengths: {len(u)} and {len(v)}"
        )
    return list(zip(u, v))


# -- Norms: the size of one vector --------------------------------------------


def p_norm(v: Vector, p: float) -> float:
    """The general p-norm of `v`: (sum of |x| ** p) ** (1 / p).

    `p` must be at least 1. Below 1 the shape stops obeying the triangle
    inequality and the result is no longer a norm at all, which is why this
    refuses rather than returning a plausible number.

    `p = math.inf` is the limit as p grows: the largest single absolute
    component, and it is computed as that limit rather than by arithmetic,
    because `x ** math.inf` overflows.

        p_norm((3, 4), 1)        -> 7.0
        p_norm((3, 4), 2)        -> 5.0
        p_norm((3, 4), math.inf) -> 4.0
    """
    if p < 1:
        raise ValueError(f"p must be at least 1 to be a norm; got {p}")
    if math.isinf(p):
        return max((abs(x) for x in v), default=0.0)
    return sum(abs(x) ** p for x in v) ** (1.0 / p)


def l1_norm(v: Vector) -> float:
    """The L1 norm: the sum of absolute values. Also called the taxicab norm."""
    return sum(abs(x) for x in v)


def l2_norm(v: Vector) -> float:
    """The L2 norm: the square root of the sum of squares.

    This is the one ordinary geometry gives you, and the only p for which the
    unit ball is round.
    """
    return math.sqrt(sum(x * x for x in v))


def linf_norm(v: Vector) -> float:
    """The L-infinity norm: the largest single absolute component."""
    return max((abs(x) for x in v), default=0.0)


# -- Distances: how far apart two vectors are ---------------------------------


def minkowski_distance(u: Vector, v: Vector, p: float) -> float:
    """The p-norm of the difference. Every distance below is a special case."""
    return p_norm([a - b for a, b in _paired(u, v)], p)


def l1_distance(u: Vector, v: Vector) -> float:
    """Manhattan distance: total disagreement, summed across the features.

    Every unit of difference costs the same wherever it happens, so ten
    features each one out costs exactly what one feature ten out costs.
    """
    return sum(abs(a - b) for a, b in _paired(u, v))


def l2_distance(u: Vector, v: Vector) -> float:
    """Euclidean distance: straight-line separation.

    Squaring makes one large disagreement cost far more than several small
    ones adding to the same total.
    """
    return math.sqrt(sum((a - b) ** 2 for a, b in _paired(u, v)))


def linf_distance(u: Vector, v: Vector) -> float:
    """Chebyshev distance: the single worst feature decides, alone.

    Every other feature is ignored entirely. That sounds like a weakness until
    you meet a tolerance check, where a part is out of specification if ANY
    dimension is out, or a two-axis machine, where the slower axis sets the
    time.
    """
    return max((abs(a - b) for a, b in _paired(u, v)), default=0.0)


# -- Angle: a similarity, not a distance --------------------------------------


def dot(u: Vector, v: Vector) -> float:
    """Day 103's dot product, repeated here so this module stands alone."""
    return sum(a * b for a, b in _paired(u, v))


def cosine_similarity(u: Vector, v: Vector) -> float:
    """The cosine of the angle between two vectors: 1 identical, 0 orthogonal.

    Day 103 derived this. It is repeated rather than re-taught. Length is
    divided out, which is the whole point: a document three times as long with
    the same word mix scores exactly 1.0.

    Undefined for the zero vector, which has no direction, so this raises
    rather than returning 0.0 and letting a silent wrong answer propagate.
    """
    nu, nv = l2_norm(u), l2_norm(v)
    if nu <= TOL or nv <= TOL:
        raise ValueError("cosine similarity is undefined for a zero vector")
    return dot(u, v) / (nu * nv)


def cosine_distance(u: Vector, v: Vector) -> float:
    """1 minus the cosine similarity.

    Widely used, useful, and NOT a metric: it fails the triangle inequality,
    which Day 103 proved and `03_metrics_and_non_metrics.py` demonstrates again
    with a concrete counter-example. Call it a dissimilarity if you want to be
    precise.
    """
    return 1.0 - cosine_similarity(u, v)


# -- Categorical and set data -------------------------------------------------


def hamming_distance(a: Sequence, b: Sequence) -> int:
    """How many positions differ. The right answer for categorical features.

    Nothing is subtracted, so the values need not be numbers: 'red' against
    'blue' is a difference of 1, exactly like 'red' against 'green'. There is
    no sense in which red is nearer to blue than to green, and any measure that
    invents one has invented data.
    """
    return sum(1 for x, y in _paired(a, b) if x != y)


def normalised_hamming(a: Sequence, b: Sequence) -> float:
    """Hamming distance as a fraction of the fields, so lengths compare."""
    if not a:
        raise ValueError("normalised Hamming needs at least one field")
    return hamming_distance(a, b) / len(a)


def jaccard_similarity(a: Iterable[Hashable], b: Iterable[Hashable]) -> float:
    """|intersection| / |union| for two sets: 1 identical, 0 disjoint.

    Two empty sets are defined here as identical, similarity 1.0. That is a
    convention rather than a derivation, and it is stated rather than hidden.
    """
    sa, sb = set(a), set(b)
    union = sa | sb
    if not union:
        return 1.0
    return len(sa & sb) / len(union)


def jaccard_distance(a: Iterable[Hashable], b: Iterable[Hashable]) -> float:
    """1 minus Jaccard similarity. Unlike cosine distance, this IS a metric."""
    return 1.0 - jaccard_similarity(a, b)


def vocabulary(*collections: Iterable[Hashable]) -> list[Hashable]:
    """The sorted union of several collections: a fixed axis order.

    Sorted rather than in encounter order, so that the binary vectors below are
    the same on every run and on every machine. A set has no order, and a
    vector built from one without sorting is a different vector each time the
    interpreter starts.
    """
    seen: set[Hashable] = set()
    for collection in collections:
        seen |= set(collection)
    return sorted(seen)


def to_binary_vector(items: Iterable[Hashable],
                     axes: Sequence[Hashable]) -> list[float]:
    """Turn a set into a 1/0 vector over a fixed list of axes.

    This is how a set gets handed to a measure that expects numbers, and it is
    where the Jaccard-against-cosine comparison becomes possible: the same two
    sets, scored both ways, on exactly the same data.
    """
    present = set(items)
    return [1.0 if axis in present else 0.0 for axis in axes]


# -- Small matrix arithmetic, in pure Python ----------------------------------


def transpose(m: Matrix) -> Matrix:
    """Rows become columns."""
    return [list(col) for col in zip(*m)]


def matmul(a: Matrix, b: Matrix) -> Matrix:
    """Day 101's matrix product: row of `a` dotted with column of `b`."""
    if len(a[0]) != len(b):
        raise DimensionMismatch(
            f"cannot multiply {len(a)}x{len(a[0])} by {len(b)}x{len(b[0])}"
        )
    bt = transpose(b)
    return [[sum(x * y for x, y in zip(row, col)) for col in bt] for row in a]


def mat_vec(m: Matrix, v: Vector) -> list[float]:
    """Matrix times column vector."""
    if len(m[0]) != len(v):
        raise DimensionMismatch(
            f"cannot apply {len(m)}x{len(m[0])} matrix to a vector of {len(v)}"
        )
    return [sum(x * y for x, y in zip(row, v)) for row in m]


def inverse(m: Matrix) -> Matrix:
    """Invert a square matrix by Gauss-Jordan elimination with partial pivoting.

    Written out rather than imported so the Mahalanobis distance below owes
    nothing to NumPy. `test_reference.py` checks it against `numpy.linalg.inv`.

    Partial pivoting -- always taking the largest available pivot -- is not
    tidiness. Without it, a small pivot divides the rest of the row by
    something near zero and the error in every later step is multiplied by
    however small it was.
    """
    n = len(m)
    if any(len(row) != n for row in m):
        raise DimensionMismatch("only a square matrix can be inverted")
    # Work on [m | I] and reduce the left half to the identity.
    aug = [list(map(float, row)) + [1.0 if i == j else 0.0 for j in range(n)]
           for i, row in enumerate(m)]
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(aug[r][col]))
        if abs(aug[pivot][col]) <= TOL:
            raise ValueError("matrix is singular: it has no inverse")
        aug[col], aug[pivot] = aug[pivot], aug[col]
        scale = aug[col][col]
        aug[col] = [x / scale for x in aug[col]]
        for row in range(n):
            if row == col:
                continue
            factor = aug[row][col]
            if factor == 0.0:
                continue
            aug[row] = [x - factor * y for x, y in zip(aug[row], aug[col])]
    return [row[n:] for row in aug]


# -- Statistics of a table of rows --------------------------------------------


def column_means(rows: Sequence[Vector]) -> list[float]:
    """The mean of each column of a table of equal-length rows."""
    if not rows:
        raise ValueError("no rows")
    n = len(rows)
    return [sum(row[j] for row in rows) / n for j in range(len(rows[0]))]


def column_stds(rows: Sequence[Vector]) -> list[float]:
    """The POPULATION standard deviation of each column: divided by n, not n-1.

    Which divisor to use is a real decision and this lab makes it explicitly.
    `n` describes the table you have; `n - 1` estimates a wider population you
    are sampling from. Scikit-learn's StandardScaler divides by `n`, and
    `numpy.std` does too unless you pass `ddof=1`, so `n` is what this lab
    uses and what the tests check against.
    """
    means = column_means(rows)
    n = len(rows)
    return [math.sqrt(sum((row[j] - means[j]) ** 2 for row in rows) / n)
            for j in range(len(rows[0]))]


def standardise(
    rows: Sequence[Vector],
    means: Sequence[float] | None = None,
    stds: Sequence[float] | None = None,
) -> list[list[float]]:
    """Subtract the column mean and divide by the column standard deviation.

    Also called the z-score. After it, every column has mean 0 and standard
    deviation 1, so a metre and a gram contribute on the same terms.

    `means` and `stds` may be supplied so that a query is standardised with the
    SAME numbers as the catalogue it is being compared against. Standardising a
    single query against itself would give a row of zeros, which is a mistake
    with a long history in production retrieval systems.

    A column with zero spread is left at 0.0 rather than dividing by zero: it
    carries no information to scale.
    """
    mu = list(means) if means is not None else column_means(rows)
    sd = list(stds) if stds is not None else column_stds(rows)
    out = []
    for row in rows:
        out.append([0.0 if sd[j] <= TOL else (row[j] - mu[j]) / sd[j]
                    for j in range(len(row))])
    return out


def covariance_matrix(rows: Sequence[Vector]) -> Matrix:
    """The population covariance matrix of a table of rows, divided by n.

    Entry (i, j) is the average product of column i's and column j's deviations
    from their own means. The diagonal is each column's variance; the
    off-diagonal is how the two move together. Day 106's eigenvectors of this
    matrix are the directions the data actually spreads along, and Mahalanobis
    distance below is what you get by measuring in those directions.
    """
    means = column_means(rows)
    n, k = len(rows), len(rows[0])
    return [[sum((row[i] - means[i]) * (row[j] - means[j]) for row in rows) / n
             for j in range(k)] for i in range(k)]


def mahalanobis_distance(u: Vector, v: Vector, cov_inverse: Matrix) -> float:
    """Euclidean distance after accounting for how the data actually varies.

    The arithmetic is one line: take the difference, and instead of dotting it
    with itself, dot it with itself THROUGH the inverse covariance matrix.

        d = sqrt( z . (cov_inverse . z) )   where z = u - v

    Substituting the identity matrix for `cov_inverse` gives back ordinary
    Euclidean distance exactly, which is the cleanest way to see what the
    covariance is doing: it re-weights the axes so that a step of one is a step
    of one standard deviation *of the data*, and it un-tilts correlated
    features so that moving along the grain of the data is cheap and moving
    across it is expensive.

    `cov_inverse` is passed in already inverted because in practice you invert
    the covariance once and then score thousands of points against it.
    """
    z = [a - b for a, b in _paired(u, v)]
    squared = dot(z, mat_vec(cov_inverse, z))
    # A covariance matrix is positive semi-definite, so `squared` is
    # non-negative in exact arithmetic. In floating point a value that should
    # be 0 can come out as -1e-17, and math.sqrt would raise on it.
    if squared < 0.0:
        if squared < -TOL:
            raise ValueError(
                f"negative squared distance ({squared}): "
                "the matrix supplied is not a valid inverse covariance"
            )
        squared = 0.0
    return math.sqrt(squared)


# -- One ranking function, with the measure as a parameter --------------------


def rank(
    query,
    candidates: dict[str, object],
    measure: Callable,
    higher_is_better: bool = False,
) -> list[tuple[str, float]]:
    """Score every candidate against the query and sort best first.

    This is the function the day is really about. Swapping Manhattan for
    cosine is ONE argument here, and the rankings move -- which is the honest
    way to see that the choice of measure is a modelling decision and not a
    detail.

    Ties are broken by the candidate's name, so the output is deterministic and
    two runs never disagree for a reason that has nothing to do with the data.
    """
    scored = [(name, float(measure(query, value)))
              for name, value in candidates.items()]
    scored.sort(key=lambda pair: (-pair[1] if higher_is_better else pair[1],
                                  pair[0]))
    return scored


def winner(query, candidates: dict[str, object], measure: Callable,
           higher_is_better: bool = False) -> str:
    """The name at the top of `rank`."""
    return rank(query, candidates, measure, higher_is_better)[0][0]

"""Tests over the reference implementation. Run from the LAB DIRECTORY:

    .venv/bin/pytest examples -q

Every float comparison states a tolerance. `measures.TOL` is 1e-12 and is used
unless a test says otherwise and says why. Integer results -- Hamming counts,
rankings, set sizes -- are compared exactly, because they are exact.

NumPy appears here as the independent answer, never as the implementation:
`measures.py` computes with `abs`, `**`, `sum` and `math.sqrt` only.
"""

from __future__ import annotations

import itertools
import math

import numpy as np
import pytest

import catalogue
import measures
from measures import TOL

Q = catalogue.QUERY
ARTICLES = catalogue.ARTICLES


# -- 1. Norms -----------------------------------------------------------------


def test_1_01_l1_norm_is_the_sum_of_absolute_values():
    assert measures.l1_norm((3.0, -4.0, 12.0)) == 19.0


def test_1_02_l2_norm_of_3_4_12_is_exactly_13():
    # 9 + 16 + 144 = 169, a perfect square, so this is exact and `==` is safe.
    assert measures.l2_norm((3.0, -4.0, 12.0)) == 13.0


def test_1_03_linf_norm_is_the_largest_absolute_component():
    assert measures.linf_norm((3.0, -4.0, 12.0)) == 12.0


def test_1_04_p_norm_reproduces_l1_l2_and_linf():
    v = (3.0, 4.0)
    assert measures.p_norm(v, 1) == measures.l1_norm(v) == 7.0
    assert measures.p_norm(v, 2) == measures.l2_norm(v) == 5.0
    assert measures.p_norm(v, math.inf) == measures.linf_norm(v) == 4.0


def test_1_05_p_norm_is_non_increasing_in_p():
    v = (3.0, 4.0)
    values = [measures.p_norm(v, p) for p in (1, 1.5, 2, 3, 4, 8, 16, 64)]
    for earlier, later in zip(values, values[1:]):
        assert later <= earlier + TOL


def test_1_06_p_norm_never_falls_below_the_largest_component():
    v = (3.0, 4.0)
    floor = measures.linf_norm(v)
    for p in (1, 1.5, 2, 3, 10, 100, math.inf):
        assert measures.p_norm(v, p) >= floor - TOL


def test_1_07_p_norm_refuses_p_below_one_because_it_is_not_a_norm():
    with pytest.raises(ValueError):
        measures.p_norm((3.0, 4.0), 0.5)


def test_1_08_every_norm_of_the_zero_vector_is_zero():
    zero = (0.0, 0.0, 0.0)
    assert measures.l1_norm(zero) == 0.0
    assert measures.l2_norm(zero) == 0.0
    assert measures.linf_norm(zero) == 0.0
    assert measures.p_norm(zero, 3) == 0.0


@pytest.mark.parametrize("norm", [measures.l1_norm, measures.l2_norm,
                                  measures.linf_norm])
def test_1_09_absolute_homogeneity(norm):
    v = catalogue.AXIOM_VECTOR
    k = catalogue.AXIOM_SCALAR
    assert abs(norm([k * x for x in v]) - abs(k) * norm(v)) <= TOL


@pytest.mark.parametrize("norm", [measures.l1_norm, measures.l2_norm,
                                  measures.linf_norm])
def test_1_10_triangle_inequality_for_norms(norm):
    v = catalogue.AXIOM_VECTOR
    w = catalogue.TRIANGLE_TRIPLE[1]
    assert norm([a + b for a, b in zip(v, w)]) <= norm(v) + norm(w) + TOL


def test_1_11_squared_euclidean_fails_absolute_homogeneity():
    """The reason 'squared distance' is not a distance, measured."""
    v = catalogue.AXIOM_VECTOR
    squared = sum(x * x for x in v)
    doubled = sum((2 * x) ** 2 for x in v)
    assert abs(doubled - 4 * squared) <= TOL
    assert abs(doubled - 2 * squared) > 1.0


def test_1_12_squared_euclidean_still_ranks_identically_to_l2():
    by_l2 = measures.rank(Q, ARTICLES, measures.l2_distance)
    by_sq = measures.rank(
        Q, ARTICLES, lambda a, b: sum((x - y) ** 2 for x, y in zip(a, b)))
    assert [n for n, _ in by_l2] == [n for n, _ in by_sq]


# -- 2. Agreement with NumPy --------------------------------------------------


@pytest.mark.parametrize("p", [1, 1.5, 2, 3, 8, math.inf])
def test_2_01_p_norm_matches_numpy_linalg_norm_ord(p):
    for v in ((3.0, 4.0), catalogue.AXIOM_VECTOR, (0.5, -0.25, 7.0, -3.5)):
        mine = measures.p_norm(v, p)
        theirs = float(np.linalg.norm(np.asarray(v), ord=p))
        assert abs(mine - theirs) <= TOL


@pytest.mark.parametrize("p, ord_", [(1, 1), (2, 2), (math.inf, np.inf)])
def test_2_02_distances_match_numpy_on_the_articles(p, ord_):
    for vec in ARTICLES.values():
        mine = measures.minkowski_distance(Q, vec, p)
        d = np.asarray(Q, dtype=float) - np.asarray(vec, dtype=float)
        assert abs(mine - float(np.linalg.norm(d, ord=ord_))) <= TOL


def test_2_03_l2_distance_matches_math_dist():
    for vec in ARTICLES.values():
        assert abs(measures.l2_distance(Q, vec)
                   - math.dist(Q, vec)) <= TOL


def test_2_04_cosine_similarity_matches_a_numpy_computation():
    for vec in ARTICLES.values():
        a = np.asarray(Q, dtype=float)
        b = np.asarray(vec, dtype=float)
        theirs = float(a @ b / (np.linalg.norm(a) * np.linalg.norm(b)))
        assert abs(measures.cosine_similarity(Q, vec) - theirs) <= TOL


# -- 3. The opening disagreement ----------------------------------------------


def test_3_01_cartogram_is_exactly_three_times_the_query():
    assert ARTICLES["Cartogram"] == tuple(3 * c for c in Q)


def test_3_02_l1_picks_aisle():
    assert measures.winner(Q, ARTICLES, measures.l1_distance) == "Aisle"


def test_3_03_l2_picks_beacon():
    assert measures.winner(Q, ARTICLES, measures.l2_distance) == "Beacon"


def test_3_04_linf_also_picks_beacon_for_a_different_reason():
    assert measures.winner(Q, ARTICLES, measures.linf_distance) == "Beacon"
    assert measures.linf_distance(Q, ARTICLES["Aisle"]) == 5.0
    assert measures.linf_distance(Q, ARTICLES["Beacon"]) == 2.0


def test_3_05_cosine_picks_cartogram():
    assert measures.winner(Q, ARTICLES, measures.cosine_similarity,
                           higher_is_better=True) == "Cartogram"


def test_3_06_three_measures_name_three_different_winners():
    picks = {
        measures.winner(Q, ARTICLES, measures.l1_distance),
        measures.winner(Q, ARTICLES, measures.l2_distance),
        measures.winner(Q, ARTICLES, measures.cosine_similarity,
                        higher_is_better=True),
    }
    assert len(picks) == 3


def test_3_07_the_exact_distances():
    assert measures.l1_distance(Q, ARTICLES["Aisle"]) == 5.0
    assert measures.l1_distance(Q, ARTICLES["Beacon"]) == 6.0
    assert measures.l1_distance(Q, ARTICLES["Cartogram"]) == 20.0
    assert measures.l2_distance(Q, ARTICLES["Aisle"]) == 5.0
    assert abs(measures.l2_distance(Q, ARTICLES["Beacon"])
               - math.sqrt(12.0)) <= TOL


def test_3_08_cosine_of_cartogram_is_exactly_one_to_tolerance():
    assert abs(measures.cosine_similarity(Q, ARTICLES["Cartogram"])
               - 1.0) <= TOL


def test_3_09_cartogram_is_the_worst_answer_under_both_l1_and_l2():
    for measure in (measures.l1_distance, measures.l2_distance):
        order = measures.rank(Q, ARTICLES, measure)
        assert order[-1][0] == "Cartogram"


def test_3_10_higher_is_better_actually_reverses_the_order():
    ascending = measures.rank(Q, ARTICLES, measures.cosine_similarity)
    descending = measures.rank(Q, ARTICLES, measures.cosine_similarity,
                               higher_is_better=True)
    assert [n for n, _ in ascending] == [n for n, _ in descending][::-1]


# -- 4. Metrics and non-metrics -----------------------------------------------


@pytest.mark.parametrize("d", [measures.l1_distance, measures.l2_distance,
                               measures.linf_distance])
def test_4_01_symmetry(d):
    x, y, _ = catalogue.TRIANGLE_TRIPLE
    assert abs(d(x, y) - d(y, x)) <= TOL


@pytest.mark.parametrize("d", [measures.l1_distance, measures.l2_distance,
                               measures.linf_distance])
def test_4_02_zero_only_when_identical(d):
    x, y, _ = catalogue.TRIANGLE_TRIPLE
    assert d(x, x) == 0.0
    assert d(x, y) > 0.0


@pytest.mark.parametrize("d", [measures.l1_distance, measures.l2_distance,
                               measures.linf_distance])
def test_4_03_triangle_inequality_in_all_six_orderings(d):
    for a, b, c in itertools.permutations(catalogue.TRIANGLE_TRIPLE):
        assert d(a, b) + d(b, c) >= d(a, c) - TOL


def test_4_04_cosine_distance_violates_the_triangle_inequality():
    east, diagonal, north = (catalogue.EAST, catalogue.DIAGONAL,
                             catalogue.NORTH)
    detour = (measures.cosine_distance(east, diagonal)
              + measures.cosine_distance(diagonal, north))
    direct = measures.cosine_distance(east, north)
    assert detour < direct - TOL
    assert abs(direct - 1.0) <= TOL
    assert abs(detour - (2.0 - math.sqrt(2.0))) <= TOL


def test_4_05_cosine_distance_is_zero_between_different_vectors():
    assert abs(measures.cosine_distance((1.0, 0.0), (2.0, 0.0))) <= TOL


def test_4_06_angular_distance_repairs_the_triangle_inequality():
    def angular(u, v):
        c = max(-1.0, min(1.0, measures.cosine_similarity(u, v)))
        return math.acos(c) / math.pi

    east, diagonal, north = (catalogue.EAST, catalogue.DIAGONAL,
                             catalogue.NORTH)
    assert (angular(east, diagonal) + angular(diagonal, north)
            >= angular(east, north) - TOL)


def test_4_07_jaccard_distance_is_a_metric_on_every_triple():
    universe = ("a", "b", "c", "d")
    subsets = [frozenset(c) for r in range(len(universe) + 1)
               for c in itertools.combinations(universe, r)]
    assert len(subsets) == 16
    for a, b, c in itertools.product(subsets, repeat=3):
        assert (measures.jaccard_distance(a, b)
                + measures.jaccard_distance(b, c)
                >= measures.jaccard_distance(a, c) - TOL)


def test_4_08_hamming_distance_is_a_metric_on_every_triple():
    strings = list(itertools.product((0, 1), repeat=4))
    for a, b, c in itertools.product(strings, repeat=3):
        assert (measures.hamming_distance(a, b)
                + measures.hamming_distance(b, c)
                >= measures.hamming_distance(a, c))


def test_4_09_cosine_distance_fails_on_many_binary_triples():
    vectors = [v for v in itertools.product((0, 1), repeat=4) if any(v)]
    violations = sum(
        1 for a, b, c in itertools.product(vectors, repeat=3)
        if (measures.cosine_distance(a, b) + measures.cosine_distance(b, c)
            < measures.cosine_distance(a, c) - TOL))
    assert violations > 0
    assert violations == 326  # measured, and asserted so a change is noticed


def test_4_10_normalised_vectors_make_cosine_and_l2_agree():
    for u, v in itertools.combinations(
            (catalogue.EAST, catalogue.DIAGONAL, catalogue.NORTH), 2):
        un = [c / measures.l2_norm(u) for c in u]
        vn = [c / measures.l2_norm(v) for c in v]
        lhs = measures.l2_distance(un, vn) ** 2
        rhs = 2.0 - 2.0 * measures.cosine_similarity(un, vn)
        assert abs(lhs - rhs) <= 1e-9


def test_4_11_cosine_similarity_refuses_the_zero_vector():
    with pytest.raises(ValueError):
        measures.cosine_similarity((0.0, 0.0), (1.0, 1.0))


# -- 5. Manhattan, Chebyshev and the shape of the question --------------------


def test_5_01_the_warehouse_displacement_gives_14_10_and_8():
    a, b = catalogue.FLOOR_FROM, catalogue.FLOOR_TO
    assert measures.l1_distance(a, b) == 14.0
    assert measures.l2_distance(a, b) == 10.0
    assert measures.linf_distance(a, b) == 8.0


def test_5_02_linf_never_exceeds_l2_never_exceeds_l1():
    rng = np.random.default_rng(107)
    for _ in range(500):
        u = rng.normal(size=5)
        v = rng.normal(size=5)
        assert (measures.linf_distance(u, v)
                <= measures.l2_distance(u, v) + TOL
                <= measures.l1_distance(u, v) + 2 * TOL)


def test_5_03_chebyshev_accepts_the_part_l1_and_l2_would_reject():
    nominal = catalogue.NOMINAL_PART
    tol_mm = catalogue.PART_TOLERANCE_MM
    a = catalogue.MEASURED_PARTS["batch-A"]
    b = catalogue.MEASURED_PARTS["batch-B"]
    assert measures.linf_distance(a, nominal) <= tol_mm + TOL
    assert measures.linf_distance(b, nominal) > tol_mm
    # ... while both L1 and L2 rank batch-B as the better part.
    assert measures.l1_distance(b, nominal) < measures.l1_distance(a, nominal)
    assert measures.l2_distance(b, nominal) > measures.l2_distance(a, nominal)


def test_5_04_l1_and_l2_disagree_about_the_two_batches():
    """The pair that shows L1 and L2 are genuinely different questions."""
    nominal = catalogue.NOMINAL_PART
    a = catalogue.MEASURED_PARTS["batch-A"]
    b = catalogue.MEASURED_PARTS["batch-B"]
    assert measures.l1_distance(a, nominal) > measures.l1_distance(b, nominal)
    assert measures.l2_distance(a, nominal) < measures.l2_distance(b, nominal)


# -- 6. Categorical and set data ----------------------------------------------


def test_6_01_hamming_counts_the_fields_that_differ():
    ref = catalogue.REFERENCE_RECORD
    got = {n: measures.hamming_distance(ref, r)
           for n, r in catalogue.CANDIDATE_RECORDS.items()}
    assert got == {"part-71": 1, "part-72": 3, "part-73": 6}


def test_6_02_normalised_hamming_is_a_fraction_of_the_fields():
    ref = catalogue.REFERENCE_RECORD
    assert measures.normalised_hamming(
        ref, catalogue.CANDIDATE_RECORDS["part-73"]) == 1.0
    assert abs(measures.normalised_hamming(
        ref, catalogue.CANDIDATE_RECORDS["part-71"]) - 1 / 6) <= TOL


def test_6_03_hamming_needs_no_arithmetic_on_the_values():
    assert measures.hamming_distance(("red", "blue"), ("green", "blue")) == 1
    assert measures.hamming_distance((None, 3, "x"), (None, 4, "x")) == 1


def test_6_04_hamming_refuses_different_lengths():
    with pytest.raises(measures.DimensionMismatch):
        measures.hamming_distance((1, 2, 3), (1, 2))


def test_6_05_on_bits_hamming_equals_l1_and_squared_l2():
    a, b = catalogue.FLAGS_A, catalogue.FLAGS_B
    h = measures.hamming_distance(a, b)
    assert h == 2
    assert abs(measures.l1_distance(a, b) - h) <= TOL
    assert abs(measures.l2_distance(a, b) ** 2 - h) <= 1e-9


def test_6_06_jaccard_and_cosine_rank_the_recipes_differently():
    q = catalogue.RECIPE_QUERY
    axes = measures.vocabulary(q, *catalogue.RECIPES.values())
    qv = measures.to_binary_vector(q, axes)
    jac = {n: measures.jaccard_similarity(q, s)
           for n, s in catalogue.RECIPES.items()}
    cos = {n: measures.cosine_similarity(qv, measures.to_binary_vector(s, axes))
           for n, s in catalogue.RECIPES.items()}
    assert max(jac, key=jac.get) == "Shortbread"
    assert max(cos, key=cos.get) == "Sachertorte"


def test_6_07_the_exact_jaccard_and_cosine_values():
    q = catalogue.RECIPE_QUERY
    axes = measures.vocabulary(q, *catalogue.RECIPES.values())
    qv = measures.to_binary_vector(q, axes)
    sach = catalogue.RECIPES["Sachertorte"]
    short = catalogue.RECIPES["Shortbread"]
    assert abs(measures.jaccard_similarity(q, sach) - 4 / 11) <= TOL
    assert abs(measures.jaccard_similarity(q, short) - 2 / 5) <= TOL
    assert abs(measures.cosine_similarity(
        qv, measures.to_binary_vector(sach, axes)) - 4 / math.sqrt(44)) <= TOL
    assert abs(measures.cosine_similarity(
        qv, measures.to_binary_vector(short, axes)) - 2 / math.sqrt(12)) <= TOL


def test_6_08_cosine_on_binary_data_is_never_below_jaccard():
    universe = tuple("abcdef")
    subsets = [frozenset(c) for r in range(1, len(universe) + 1)
               for c in itertools.combinations(universe, r)]
    axes = list(universe)
    for a, b in itertools.combinations(subsets, 2):
        jac = measures.jaccard_similarity(a, b)
        cos = measures.cosine_similarity(measures.to_binary_vector(a, axes),
                                         measures.to_binary_vector(b, axes))
        assert cos >= jac - TOL


def test_6_09_jaccard_of_two_empty_sets_is_the_documented_convention():
    assert measures.jaccard_similarity(set(), set()) == 1.0
    assert measures.jaccard_distance(set(), set()) == 0.0


def test_6_10_jaccard_of_disjoint_sets_is_zero():
    assert measures.jaccard_similarity({1, 2}, {3, 4}) == 0.0


def test_6_11_vocabulary_is_sorted_so_the_axes_are_stable():
    axes = measures.vocabulary({"pear", "apple"}, {"fig"})
    assert axes == ["apple", "fig", "pear"]


def test_6_12_to_binary_vector_marks_exactly_the_members():
    axes = ["a", "b", "c"]
    assert measures.to_binary_vector({"a", "c"}, axes) == [1.0, 0.0, 1.0]


# -- 7. Matrices, covariance and Mahalanobis ---------------------------------


def test_7_01_covariance_of_the_readings_is_exactly_7_5_and_7():
    assert measures.covariance_matrix(catalogue.SENSOR_READINGS) == [
        [7.5, 7.0], [7.0, 7.5]]


def test_7_02_covariance_matches_numpy_with_bias_true():
    mine = np.asarray(measures.covariance_matrix(catalogue.SENSOR_READINGS))
    theirs = np.cov(np.asarray(catalogue.SENSOR_READINGS, dtype=float),
                    rowvar=False, bias=True)
    assert np.allclose(mine, theirs, atol=TOL)


def test_7_03_inverse_matches_numpy_linalg_inv():
    for m in ([[7.5, 7.0], [7.0, 7.5]],
              [[2.0, 0.0, 1.0], [1.0, 3.0, 2.0], [1.0, 0.0, 4.0]],
              [[1.0, 0.0], [0.0, 1.0]]):
        mine = np.asarray(measures.inverse(m))
        assert np.allclose(mine, np.linalg.inv(np.asarray(m)), atol=1e-10)


def test_7_04_inverse_times_original_is_the_identity():
    m = [[7.5, 7.0], [7.0, 7.5]]
    product = measures.matmul(m, measures.inverse(m))
    for i, row in enumerate(product):
        for j, value in enumerate(row):
            assert abs(value - (1.0 if i == j else 0.0)) <= TOL


def test_7_05_a_singular_matrix_refuses_to_invert():
    with pytest.raises(ValueError):
        measures.inverse([[1.0, 2.0], [2.0, 4.0]])


def test_7_06_the_two_probes_are_the_same_euclidean_distance_from_the_mean():
    mean = measures.column_means(catalogue.SENSOR_READINGS)
    a = measures.l2_distance(catalogue.PROBE_ALONG, mean)
    b = measures.l2_distance(catalogue.PROBE_ACROSS, mean)
    assert abs(a - b) <= TOL
    assert abs(a - math.sqrt(18.0)) <= TOL


def test_7_07_mahalanobis_tells_them_apart():
    mean = measures.column_means(catalogue.SENSOR_READINGS)
    inv = measures.inverse(measures.covariance_matrix(
        catalogue.SENSOR_READINGS))
    along = measures.mahalanobis_distance(catalogue.PROBE_ALONG, mean, inv)
    across = measures.mahalanobis_distance(catalogue.PROBE_ACROSS, mean, inv)
    assert abs(across - 6.0) <= TOL
    assert abs(along - math.sqrt(9.0 / 7.25)) <= TOL
    assert across / along > 5.0


def test_7_08_the_same_two_probes_through_numpy_agree_within_tolerance():
    """Two correct inverses disagree in the last bit. The tolerance earns
    its keep here: `== 6.0` passes for one route and fails for the other."""
    data = np.asarray(catalogue.SENSOR_READINGS, dtype=float)
    inv_np = np.linalg.inv(np.cov(data, rowvar=False, bias=True))
    z = np.asarray(catalogue.PROBE_ACROSS)
    theirs = float(math.sqrt(z @ inv_np @ z))
    mine = measures.mahalanobis_distance(
        catalogue.PROBE_ACROSS,
        measures.column_means(catalogue.SENSOR_READINGS),
        measures.inverse(measures.covariance_matrix(
            catalogue.SENSOR_READINGS)))
    assert abs(mine - theirs) <= TOL
    assert abs(theirs - 6.0) <= TOL


def test_7_09_mahalanobis_with_the_identity_is_euclidean():
    identity = [[1.0, 0.0], [0.0, 1.0]]
    for probe in ((3.0, 3.0), (3.0, -3.0), (-2.5, 4.75), (0.0, 0.0)):
        assert abs(measures.mahalanobis_distance(probe, (0.0, 0.0), identity)
                   - measures.l2_distance(probe, (0.0, 0.0))) <= TOL


def test_7_10_mahalanobis_agrees_with_the_eigen_decomposition():
    """Day 106's eigenvectors are the axes Mahalanobis measures along."""
    cov = measures.covariance_matrix(catalogue.SENSOR_READINGS)
    values, vectors = np.linalg.eigh(np.asarray(cov))
    inv = measures.inverse(cov)
    for probe in (catalogue.PROBE_ALONG, catalogue.PROBE_ACROSS,
                  (1.0, -4.0), (2.5, 0.25)):
        z = np.asarray(probe)
        by_hand = math.sqrt(sum(
            float(z @ vector) ** 2 / value
            for value, vector in zip(values, vectors.T)))
        assert abs(by_hand - measures.mahalanobis_distance(
            probe, (0.0, 0.0), inv)) <= 1e-9


def test_7_11_the_eigenvalues_are_half_and_fourteen_and_a_half():
    values = sorted(np.linalg.eigvalsh(
        np.asarray(measures.covariance_matrix(catalogue.SENSOR_READINGS))))
    assert abs(values[0] - 0.5) <= 1e-12
    assert abs(values[1] - 14.5) <= 1e-12


def test_7_12_mahalanobis_is_symmetric():
    inv = measures.inverse(measures.covariance_matrix(
        catalogue.SENSOR_READINGS))
    a, b = catalogue.PROBE_ALONG, catalogue.PROBE_ACROSS
    assert abs(measures.mahalanobis_distance(a, b, inv)
               - measures.mahalanobis_distance(b, a, inv)) <= TOL


# -- 8. Standardising ---------------------------------------------------------


def test_8_01_standardised_columns_have_mean_zero_and_sd_one():
    rows = list(catalogue.BEARINGS.values())
    z = measures.standardise(rows)
    for mean in measures.column_means(z):
        assert abs(mean) <= 1e-12
    for sd in measures.column_stds(z):
        assert abs(sd - 1.0) <= 1e-12


def test_8_02_column_stds_use_the_population_divisor_n():
    rows = list(catalogue.BEARINGS.values())
    mine = measures.column_stds(rows)
    theirs = np.asarray(rows, dtype=float).std(axis=0)
    assert np.allclose(np.asarray(mine), theirs, atol=TOL)
    # ... and NOT the n-1 divisor, which is a visibly different number.
    sample = np.asarray(rows, dtype=float).std(axis=0, ddof=1)
    assert not np.allclose(np.asarray(mine), sample, atol=1e-6)


def test_8_03_raw_euclidean_picks_the_unusable_bearing():
    assert measures.winner(catalogue.BEARING_QUERY, catalogue.BEARINGS,
                           measures.l2_distance) == "R"


def test_8_04_standardising_changes_the_winner():
    rows = list(catalogue.BEARINGS.values())
    means = measures.column_means(rows)
    stds = measures.column_stds(rows)
    q = measures.standardise([catalogue.BEARING_QUERY], means, stds)[0]
    scaled = {n: measures.standardise([v], means, stds)[0]
              for n, v in catalogue.BEARINGS.items()}
    assert measures.winner(q, scaled, measures.l2_distance) == "P"


def test_8_05_the_bore_column_contributes_almost_nothing_before_scaling():
    q = catalogue.BEARING_QUERY
    for row in catalogue.BEARINGS.values():
        bore = (q[0] - row[0]) ** 2
        mass = (q[1] - row[1]) ** 2
        if bore + mass > 0:
            assert bore / (bore + mass) < 1e-4


def test_8_06_a_unit_change_alone_flips_the_ranking():
    q, cands = catalogue.BEARING_QUERY, catalogue.BEARINGS
    metres = measures.winner(q, cands, measures.l2_distance)
    micro_q = (q[0] * 1e6, q[1])
    micro = {n: (v[0] * 1e6, v[1]) for n, v in cands.items()}
    assert metres == "R"
    assert measures.winner(micro_q, micro, measures.l2_distance) == "P"


def test_8_07_standardise_leaves_a_constant_column_at_zero():
    rows = [(1.0, 5.0), (2.0, 5.0), (3.0, 5.0)]
    z = measures.standardise(rows)
    assert [row[1] for row in z] == [0.0, 0.0, 0.0]


def test_8_08_standardising_a_query_against_itself_would_give_zeros():
    """The mistake the API's `means`/`stds` arguments exist to prevent."""
    alone = measures.standardise([catalogue.BEARING_QUERY])[0]
    assert alone == [0.0, 0.0]


def test_8_09_cosine_is_not_invariant_to_a_column_unit_change():
    q, cands = catalogue.BEARING_QUERY, catalogue.BEARINGS
    a = measures.rank(q, cands, measures.cosine_similarity,
                      higher_is_better=True)
    micro_q = (q[0] * 1e6, q[1])
    micro = {n: (v[0] * 1e6, v[1]) for n, v in cands.items()}
    b = measures.rank(micro_q, micro, measures.cosine_similarity,
                      higher_is_better=True)
    assert [n for n, _ in a] != [n for n, _ in b]


def test_8_10_cosine_is_invariant_to_scaling_a_whole_vector():
    q, cands = catalogue.BEARING_QUERY, catalogue.BEARINGS
    a = measures.rank(q, cands, measures.cosine_similarity,
                      higher_is_better=True)
    doubled = {n: tuple(7.5 * c for c in v) for n, v in cands.items()}
    b = measures.rank(q, doubled, measures.cosine_similarity,
                      higher_is_better=True)
    assert [n for n, _ in a] == [n for n, _ in b]
    for (_, x), (_, y) in zip(a, b):
        assert abs(x - y) <= TOL


# -- 9. The ranking function and the guard rails ------------------------------


def test_9_01_rank_returns_every_candidate_once():
    order = measures.rank(Q, ARTICLES, measures.l1_distance)
    assert sorted(n for n, _ in order) == sorted(ARTICLES)


def test_9_02_rank_is_sorted_ascending_for_a_distance():
    scores = [s for _, s in measures.rank(Q, ARTICLES, measures.l2_distance)]
    assert scores == sorted(scores)


def test_9_03_rank_is_sorted_descending_for_a_similarity():
    scores = [s for _, s in measures.rank(Q, ARTICLES,
                                          measures.cosine_similarity,
                                          higher_is_better=True)]
    assert scores == sorted(scores, reverse=True)


def test_9_04_ties_break_by_name_so_the_output_is_deterministic():
    candidates = {"zulu": (1.0, 1.0), "alpha": (1.0, 1.0),
                  "mike": (1.0, 1.0)}
    order = [n for n, _ in measures.rank((0.0, 0.0), candidates,
                                         measures.l2_distance)]
    assert order == ["alpha", "mike", "zulu"]


def test_9_05_swapping_the_measure_is_one_argument():
    """The claim the lab is built on, stated as a test."""
    picks = {name: measures.winner(Q, ARTICLES, measure, higher)
             for name, (measure, higher) in {
                 "l1": (measures.l1_distance, False),
                 "l2": (measures.l2_distance, False),
                 "cos": (measures.cosine_similarity, True),
             }.items()}
    assert picks == {"l1": "Aisle", "l2": "Beacon", "cos": "Cartogram"}


@pytest.mark.parametrize("fn", [measures.l1_distance, measures.l2_distance,
                                measures.linf_distance, measures.dot])
def test_9_06_comparing_different_lengths_raises(fn):
    with pytest.raises(measures.DimensionMismatch):
        fn((1.0, 2.0, 3.0), (1.0, 2.0))


def test_9_07_dimension_mismatch_is_catchable_as_a_value_error():
    with pytest.raises(ValueError):
        measures.l2_distance((1.0,), (1.0, 2.0))


def test_9_08_mat_vec_and_matmul_agree_with_numpy():
    m = [[2.0, -1.0, 0.5], [0.0, 3.0, 1.0]]
    v = [1.0, 2.0, -4.0]
    assert np.allclose(np.asarray(measures.mat_vec(m, v)),
                       np.asarray(m) @ np.asarray(v), atol=TOL)
    n = [[1.0, 0.0], [2.0, -1.0], [0.5, 4.0]]
    assert np.allclose(np.asarray(measures.matmul(m, n)),
                       np.asarray(m) @ np.asarray(n), atol=TOL)


def test_9_09_transpose_swaps_rows_and_columns():
    assert measures.transpose([[1, 2, 3], [4, 5, 6]]) == [[1, 4], [2, 5],
                                                          [3, 6]]


def test_9_10_mahalanobis_refuses_a_matrix_that_is_not_a_covariance_inverse():
    bad = [[-1.0, 0.0], [0.0, -1.0]]
    with pytest.raises(ValueError):
        measures.mahalanobis_distance((1.0, 1.0), (0.0, 0.0), bad)

lesson_id: D107
day: 107
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-107-norms-distances-and-similarity-measures
  - 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_three_measures_three_winners.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_the_p_norm_family.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_metrics_and_non_metrics.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_choosing_by_the_shape_of_the_data.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_mahalanobis_distance.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_scaling_changes_the_answer.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 . -name '.venv' -prune -o -type d -name '__pycache__' -print -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-17'
executed_on: 'macOS 26.5.2 (Apple Silicon, arm64), Python 3.14.0, numpy 2.5.2, pytest 9.1.1, bash 3.2.57 — bash tests/run_tests.sh -> 98 checks, 0 failure(s), exit 0; pytest examples -> 105 passed; pytest starter -> 1 passed, 71 skipped on an untouched checkout, and 72 passed against a fully solved copy of starter/ kept outside the lab. All six reference scripts exit 0 with every internal assertion holding. Network is needed once to install numpy and pytest; nothing else in the lab opens a socket, and every dataset is a literal table in catalogue.py rather than a download — section 7 of the harness greps the sources for network calls, asserts that no data file exists anywhere under the lab, and also greps both copies of measures.py to confirm neither imports NumPy, which is what makes agreeing with numpy.linalg.norm evidence rather than a tautology. Section 6 re-runs the harness with one expectation deliberately swapped for the naive belief that cosine distance satisfies the triangle inequality, 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. The harness was verified to exit 0 both with the lab-local .venv present and with it absent (PYTEST pointed elsewhere); every find in section 7 prunes .venv first, because numpy ships 113 __pycache__ directories inside it and the README instructs the reader to create it. Four measured results are asserted rather than smoothed over. (1) The same Mahalanobis distance comes out as exactly 6.0 through the lab hand-written Gauss-Jordan inverse and as 5.999999999999999 through numpy.linalg.inv — two correct routes disagreeing by 8.882e-16 — and both are asserted against a tolerance of 1e-12 rather than one being declared right. (2) Cosine distance violates the triangle inequality on 326 of the 3375 triples of non-zero 4-bit vectors, while Jaccard distance and Hamming distance survive all 4096 triples of their own exhaustive sweeps. (3) Standardising the bearing catalogue moves exactly two of the six parts, and they are the two the decision is between: P and R swap first for third. (4) Mahalanobis on the RAW bearing numbers gives a THIRD ranking, U first and P second, not the standardised one — the lab brief expected it to reproduce the standardised answer and it does not, because Mahalanobis also removes the +0.7979 correlation between bore and mass while standardising only rescales each column. That disagreement is reported and explained rather than the claim being trimmed to fit. The one randomised demonstration uses a SEEDED numpy.random.default_rng(107) over 2000 catalogues; the code asserts a range (35 to 75 per cent of winners change) because NumPy does not guarantee stream stability across versions, while the harness additionally records the exact observed figure of 1090, and expected-output/FIELDS.md explains which of the two claims is the durable one.'

# Dependencies for the Day 107 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 independent answer. `numpy.linalg.norm(v, ord=p)` IS the p-norm family this lab implements by hand, `numpy.cov` is the covariance matrix, `numpy.linalg.inv` is the inverse, and `numpy.linalg.eigh` supplies Day 106's eigenvectors. It also provides the seeded generator for the one randomised demonstration. |
| `pytest` | `9.1.1` | The test runner from Days 071–074. Nothing new here except what it is pointed at. |

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

`examples/measures.py` and `starter/measures.py` compute every norm, distance,
similarity, mean, standard deviation, covariance and matrix inverse with
`abs`, `**`, `sum`, `max` and `math.sqrt`. Nothing in either file imports
NumPy, and section 7 of `tests/run_tests.sh` greps both files to check that it
stays that way.

That is not stylistic. The lab's central evidence is that a hand-written
`p_norm` agrees with `numpy.linalg.norm(v, ord=p)` to within 1e-12 across six
values of `p`, and that a hand-written Gauss-Jordan inverse agrees with
`numpy.linalg.inv`. If the hand-written version were built out of NumPy calls,
both comparisons would be NumPy checking itself, and would prove nothing.

The one place the two genuinely disagree is preserved rather than smoothed
over: the Mahalanobis distance across the grain of the sensor data comes out as
exactly `6.0` through Gauss-Jordan and `5.999999999999999` through LAPACK.
`expected-output/FIELDS.md` explains it.

## Why the versions are pinned

They are *checked* rather than assumed. Section 1 of `tests/run_tests.sh` reads
the installed version of each package and compares it against this file, so a
mismatch is reported at the top of the run instead of surfacing later as a
confusing difference in output.

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

1. **The seeded sweep.** `numpy.random.default_rng(107)` produces a
   reproducible stream on a given NumPy build, and NumPy's own documentation
   declines to guarantee that the stream survives a version change. So the code
   asserts a *range* — the ranking winner changes in between 35 and 75 per cent
   of 2000 random catalogues — while the harness additionally records the exact
   observed figure of 1090. If only the exact figure moves, nothing is broken.

2. **Last-bit floating point.** The `5.999999999999999` above came out of this
   LAPACK build. A different one could produce `6.0`, or
   `6.000000000000001`. The lab asserts that both routes land within 1e-12 of
   6, which is the claim that actually matters, and prints both values so a
   reader can see them rather than take the claim on trust.

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

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

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

## 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.**

Every dataset is written out as a literal table in `catalogue.py`: four term
counts, six categorical records, two ingredient lists, eight sensor readings
and six bearings. That was a deliberate choice. A lab that downloads a dataset
is a lab that breaks on a train, depends on a URL serving the same bytes next
year, and hides its own test data behind a fetch so a reader cannot see what is
in it without running the code. Section 7 of `tests/run_tests.sh` greps every
file under `examples/` and `starter/` for the patterns that would indicate a
socket being opened, and also asserts that no data file exists anywhere under
the lab.

## Running without a lab-local environment

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

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

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

## What you would give up without NumPy

Less than on most days, and it is worth being precise about which half.

All seventeen functions in exercise 1 need only the standard library, and so do
all twenty-five predictions in exercise 2. Every claim the day makes about
*which measure wins* — three winners on one query, Chebyshev accepting the part
L1 rejects, Jaccard and cosine disagreeing on the same two sets, Mahalanobis
separating two Euclidean-identical points, standardising changing the ranking —
is computed entirely in `measures.py` and would still run.

What you lose is the *corroboration*: `numpy.linalg.norm(v, ord=p)` confirming
the p-norm family, `numpy.cov` confirming the covariance, `numpy.linalg.inv`
confirming the inverse, `numpy.linalg.eigh` connecting Mahalanobis back to Day
106's eigenvectors, and the seeded 2000-catalogue sweep showing the scaling
effect is not a property of six hand-picked rows. Those are the checks that
make the lab's numbers evidence rather than assertion, and the lab does not
pretend they are optional.

numpy==2.5.2
pytest==9.1.1

# Day 107 lab — Choose Your Distance on Purpose

One idea holds this whole lab together:

> **"Distance" is not one thing. It is a family, and choosing a member is a
> modelling decision with consequences you can see.**

Everything below is a consequence of that sentence. You will implement five
measures from arithmetic you can check on paper, watch three of them name three
different winners on the same four numbers, prove that one of the most popular
of them is not a metric, and then find that the answer changes again when you
divide two columns by their own standard deviations.

Work in order; each exercise uses the one before it.

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

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

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

---

## The data

Nothing is downloaded. `catalogue.py` holds five small hand-written datasets,
every number small enough to check without a calculator. Read it first — every
exercise refers to it.

| Dataset | What it is for |
| --- | --- |
| `QUERY` and `ARTICLES` | The opening disagreement: four term counts, three articles, three different winners. |
| `FLOOR_FROM` / `FLOOR_TO`, `MEASURED_PARTS` | Where Manhattan, Euclidean and Chebyshev are each the only correct answer. |
| `REFERENCE_RECORD`, `CANDIDATE_RECORDS`, `FLAGS_A/B` | Categorical and binary data, which is what Hamming is for. |
| `RECIPE_QUERY`, `RECIPES` | Sets, where Jaccard and cosine give opposite answers. |
| `SENSOR_READINGS`, `PROBE_ALONG`, `PROBE_ACROSS` | Two points Euclidean distance cannot tell apart and Mahalanobis can. |
| `BEARINGS`, `BEARING_QUERY` | Two features in mismatched units, where scaling decides the winner. |

There is one seeded random generator in the lab,
`numpy.random.default_rng(107)` in `examples/06_scaling_changes_the_answer.py`,
and it is used only to show that the scaling effect is not a property of the
six hand-picked bearings. Every number asserted to the last decimal place comes
from the literal tables.

---

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

**1. A distance shrinks as things get more alike. A similarity grows.**

`l1_distance`, `l2_distance`, `linf_distance`, `hamming_distance` and
`mahalanobis_distance` are distances: **smaller is better**.
`cosine_similarity` and `jaccard_similarity` are similarities: **larger is
better**.

Nothing in this module guesses which you meant. `rank` takes an explicit
`higher_is_better` flag, and getting it backwards returns the *worst* match at
the top of your results with no error message anywhere. This is not a
hypothetical: it is one of the most common bugs in a first retrieval system,
and it is invisible because the output still looks like a ranked list.

**2. Use only the standard library inside `measures.py`.**

`math`, `abs`, `sum`, `max`, `sorted` and set operations are everything you
need. NumPy appears in the *tests*, where it is the independent answer —
`numpy.linalg.norm(v, ord=p)` is exactly the p-norm family you are about to
write, and agreeing with it to 1e-12 means something only if your answer was
not NumPy's answer to begin with.

---

## Exercise 1 — the measures (`measures.py`)

Seventeen functions. Each has a docstring with a worked example.

| # | Function | The one thing to get right |
| --- | --- | --- |
| 1.1 | `l1_norm` | Absolute values, then sum. |
| 1.2 | `l2_norm` | `math.sqrt` of the sum of squares. |
| 1.3 | `linf_norm` | `max(..., default=0.0)` — the empty case raises without it. |
| 1.4 | `p_norm` | `p = math.inf` is the LIMIT, not arithmetic. `p < 1` must raise. |
| 1.5 | `l1_distance` | Use `_paired`, not bare `zip`. |
| 1.6 | `l2_distance` | |
| 1.7 | `linf_distance` | |
| 1.8 | `cosine_similarity` | The zero vector has no direction: raise, do not return 0. |
| 1.9 | `hamming_distance` | Count, do not subtract. Return an `int`. |
| 1.10 | `jaccard_similarity` | Two empty sets are 1.0 by convention. State it, do not divide by zero. |
| 1.11 | `to_binary_vector` | Floats, not booleans. |
| 1.12 | `column_means` | |
| 1.13 | `column_stds` | Population divisor `n`, not `n - 1`. |
| 1.14 | `standardise` | Use the supplied `means`/`stds` when given. A zero-spread column stays 0.0. |
| 1.15 | `covariance_matrix` | Divide by `n`. The answer on the sensor readings is exactly `[[7.5, 7.0], [7.0, 7.5]]`. |
| 1.16 | `mahalanobis_distance` | `sqrt(z . (cov_inverse . z))`. Clamp a tiny negative to 0. |
| 1.17 | `rank` | Sort on a TUPLE so ties break by name. |

Written for you: `_paired`, `dot`, `minkowski_distance`, `cosine_distance`,
`normalised_hamming`, `jaccard_distance`, `vocabulary`, `transpose`, `matmul`,
`mat_vec`, `inverse` and `winner`. Gauss-Jordan elimination is a day of its
own; you built matrix multiplication on Day 101.

### The order that will hurt least

1. `l1_norm`, `l2_norm`, `linf_norm` — three one-liners, and the tests for the
   distances depend on nothing else.
2. `p_norm` — then check that it reproduces all three.
3. The three distances, then `cosine_similarity`.
4. `hamming_distance`, `jaccard_similarity`, `to_binary_vector`.
5. `rank`. Everything above becomes visible the moment this works.
6. `column_means`, `column_stds`, `standardise`.
7. `covariance_matrix`, then `mahalanobis_distance`.

---

## Exercise 2 — the predictions (`answers.py`)

Twenty-five values. **Write your answer before you run anything.** A prediction
you got wrong teaches more than a result you read, and every one of these can
be worked out on paper or reasoned about from the docstrings.

Four of them are worth thinking about rather than computing:

- **2.12** — does the p-norm rise or fall as `p` grows? Try `p = 1` and
  `p = 100` on `(3, 4)` in your head before you decide.
- **2.13** — which axiom does *squared* Euclidean distance break? There are
  four candidates and only one survives contact with doubling a vector.
- **2.15** — does the triangle inequality hold for cosine distance? The triple
  is `(1, 0)`, `(1, 1)`, `(0, 1)`, and the answer is the reason Day 103's
  result matters.
- **2.22** — are the two probe points the same Euclidean distance from the
  mean? Look at the two coordinates rather than reaching for a calculator.

---

## When you are finished

```bash
.venv/bin/pytest starter -q          # should say 72 passed
bash tests/run_tests.sh              # the whole lab, 60 checks
```

Then read `examples/`, in order 01 to 06. The reference implementation is
`examples/measures.py`, and the six scripts are the argument the day makes:
three winners, the unit balls, the metric axioms, the four data shapes,
Mahalanobis, and the scaling that changes the answer.

Reading them before you have attempted the exercises is allowed and is a waste
of the exercise.

"""Exercise 2 -- twenty-five predictions. Replace each `None` with your answer.

These are not busy-work. Every one of them is a value you can work out on
paper, or a judgement the day turns on, and writing the answer down BEFORE you
run the code is the difference between reading a result and predicting one.

A `None` is SKIPPED by the test suite, not failed. Answer the ones you are
sure of, run the tests, then come back for the rest.

Numeric answers are compared with a tolerance of 1e-9 unless stated otherwise,
so `5` and `5.0` are both fine and you need not type more than four decimals
where four are asked for.

The data is in `catalogue.py`. Read it -- everything here refers to it.
"""

# -- The opening disagreement -------------------------------------------------
#
# QUERY   = (4, 3, 2, 1) over the terms (norm, distance, vector, cluster)
# Aisle     = (4, 3, 2, 6)
# Beacon    = (6, 1, 4, 1)
# Cartogram = (12, 9, 6, 3)

# 2.1  The L1 (Manhattan) distance from QUERY to Aisle.
L1_QUERY_TO_AISLE = None

# 2.2  The L1 distance from QUERY to Beacon.
L1_QUERY_TO_BEACON = None

# 2.3  The L2 (Euclidean) distance from QUERY to Aisle. It is a whole number.
L2_QUERY_TO_AISLE = None

# 2.4  The L-infinity (Chebyshev) distance from QUERY to Beacon.
LINF_QUERY_TO_BEACON = None

# 2.5  The cosine similarity between QUERY and Cartogram, to four decimals.
#      Look at the two vectors before you reach for a calculator.
COSINE_QUERY_TO_CARTOGRAM = None

# 2.6  Which article does L1 rank first? A string: "Aisle", "Beacon" or
#      "Cartogram".
L1_WINNER = None

# 2.7  Which article does L2 rank first?
L2_WINNER = None

# 2.8  Which article does cosine similarity rank first?
COSINE_WINNER = None

# -- The p-norm family --------------------------------------------------------

# 2.9  p_norm((3.0, 4.0), 1)
P_NORM_3_4_AT_1 = None

# 2.10 p_norm((3.0, 4.0), 2)
P_NORM_3_4_AT_2 = None

# 2.11 p_norm((3.0, 4.0), math.inf)
P_NORM_3_4_AT_INF = None

# 2.12 As p rises from 1 towards infinity, does the p-norm of a fixed vector
#      rise, fall, or stay the same? One of the strings "rise", "fall",
#      "stay the same".
P_NORM_AS_P_RISES = None

# 2.13 Which of the four norm axioms does SQUARED Euclidean distance break?
#      One of the strings "non-negativity", "zero only at zero",
#      "absolute homogeneity", "triangle inequality".
AXIOM_SQUARED_EUCLIDEAN_BREAKS = None

# -- Metrics ------------------------------------------------------------------

# 2.14 cosine_distance(EAST, NORTH), where EAST = (1, 0) and NORTH = (0, 1).
COSINE_DISTANCE_EAST_NORTH = None

# 2.15 Does the triangle inequality hold for cosine distance on the triple
#      EAST, DIAGONAL, NORTH? True or False.
COSINE_TRIANGLE_HOLDS = None

# 2.16 Is 1 - Jaccard similarity a metric? True or False.
JACCARD_DISTANCE_IS_A_METRIC = None

# -- Categorical and set data -------------------------------------------------

# 2.17 The Hamming distance from REFERENCE_RECORD to part-73, out of 6 fields.
HAMMING_REFERENCE_TO_PART_73 = None

# 2.18 The Hamming distance between FLAGS_A and FLAGS_B.
HAMMING_FLAGS = None

# 2.19 Which recipe does JACCARD similarity prefer? "Sachertorte" or
#      "Shortbread".
JACCARD_RECIPE_WINNER = None

# 2.20 Which recipe does COSINE similarity prefer, on the same two sets?
COSINE_RECIPE_WINNER = None

# -- Mahalanobis --------------------------------------------------------------

# 2.21 The population covariance matrix of SENSOR_READINGS, as a list of lists.
#      Every entry is a clean number; work it out rather than guessing.
COVARIANCE_OF_READINGS = None

# 2.22 Are PROBE_ALONG and PROBE_ACROSS the same Euclidean distance from the
#      mean of SENSOR_READINGS? True or False.
PROBES_EQUIDISTANT_UNDER_EUCLIDEAN = None

# 2.23 The Mahalanobis distance from the mean to PROBE_ACROSS = (3, -3).
#      It is a whole number.
MAHALANOBIS_TO_PROBE_ACROSS = None

# -- Scaling ------------------------------------------------------------------

# 2.24 Which bearing wins on the RAW numbers, bore in metres and mass in
#      grams? A single letter as a string.
RAW_BEARING_WINNER = None

# 2.25 Which bearing wins after both columns are standardised?
STANDARDISED_BEARING_WINNER = None

"""The data this lab argues over. Written by hand, on purpose.

Nothing here is downloaded and nothing here is random. Every number is small
enough to check on paper, and each dataset was chosen because it makes exactly
one measure look right and the others look wrong -- which is the point of the
day.

There IS a seeded random generator in this lab, in
`06_scaling_changes_the_answer.py`, used to show that the scaling effect is not
a property of these six hand-picked parts. It is seeded with
`numpy.random.default_rng(107)` and every claim made about it is structural (a
count, a direction, a percentage floor) rather than a specific digit, because
NumPy does not promise that a generator's exact stream survives a version
change. Everything asserted to the last decimal place in this lab comes from the
literal tables below.
"""

from __future__ import annotations

# ---------------------------------------------------------------------------
# 1. The opening disagreement: three articles, one query, three winners.
# ---------------------------------------------------------------------------
#
# A tiny help-centre search. Each article is counted over four terms, and the
# query is the reader's own short note. Raw counts, not frequencies -- which is
# exactly what makes the three measures disagree.

TERMS = ("norm", "distance", "vector", "cluster")

QUERY = (4, 3, 2, 1)

# "Aisle" mentions cluster six times where the query mentions it once, and
# matches the other three terms exactly. One big disagreement, nothing else.
#
# "Beacon" is a little off on three of the four terms and exact on the fourth.
# Three small disagreements adding to more in total than Aisle's single one.
#
# "Cartogram" is the query's profile at exactly three times the length: a long
# article on precisely this topic. Its direction is identical, so its cosine
# similarity is exactly 1.0, and its raw counts are further away than anything
# else here.
ARTICLES: dict[str, tuple[int, ...]] = {
    "Aisle": (4, 3, 2, 6),
    "Beacon": (6, 1, 4, 1),
    "Cartogram": (12, 9, 6, 3),
}

# ---------------------------------------------------------------------------
# 2. Chebyshev, and where a single worst component decides.
# ---------------------------------------------------------------------------
#
# One displacement, in metres, across a warehouse floor laid out on aisles.
# The same two points, three operationally different answers:
#
#   L1  = 14  a picker who must walk the aisles, one axis at a time
#   L2  = 10  a drone that can fly the diagonal
#   Linf = 8  a two-axis gantry whose motors run at once, so the slower axis
#             alone sets the finishing time
#
# None of the three is the "real" distance. Each is the real distance for a
# different machine.

FLOOR_FROM = (0.0, 0.0)
FLOOR_TO = (6.0, 8.0)

# A machined part and its nominal dimensions, in millimetres. The part is
# rejected if ANY dimension is out by more than the tolerance -- which is a
# Chebyshev ball, and nothing else.
NOMINAL_PART = (40.00, 25.00, 12.00, 6.00)
PART_TOLERANCE_MM = 0.05
MEASURED_PARTS: dict[str, tuple[float, ...]] = {
    # Four dimensions each a little out. Total error is large, worst is small.
    "batch-A": (40.04, 24.96, 12.04, 5.96),
    # Three dimensions perfect, one badly out. Total error is smaller.
    "batch-B": (40.00, 25.00, 12.00, 6.09),
}

# ---------------------------------------------------------------------------
# 3. Hamming, for data with no arithmetic in it.
# ---------------------------------------------------------------------------
#
# Six categorical fields from a parts register. Subtracting "steel" from
# "brass" is not a smaller number than subtracting "steel" from "nylon"; it is
# not a number at all. Hamming counts the fields that differ and refuses to
# invent an ordering.

FIELDS = ("material", "finish", "thread", "grade", "colour", "origin")

REFERENCE_RECORD = ("steel", "zinc", "M8", "8.8", "silver", "IN")
CANDIDATE_RECORDS: dict[str, tuple[str, ...]] = {
    "part-71": ("steel", "zinc", "M8", "8.8", "black", "IN"),
    "part-72": ("brass", "zinc", "M8", "10.9", "silver", "DE"),
    "part-73": ("nylon", "plain", "M6", "4.6", "white", "CN"),
}

# The same measure on bits, which is where Hamming was defined: two 8-bit
# feature flags from the same register.
FLAGS_A = (1, 0, 1, 1, 0, 0, 1, 0)
FLAGS_B = (1, 0, 0, 1, 0, 1, 1, 0)

# ---------------------------------------------------------------------------
# 4. Jaccard against cosine on the same set-like data.
# ---------------------------------------------------------------------------
#
# Ingredient lists. The query has four ingredients.
#
#   "Sachertorte" contains ALL FOUR of them, plus seven more.
#   "Shortbread"  shares two of the four and has one ingredient of its own.
#
# Cosine, which divides by the square root of the two sizes, prefers
# Sachertorte: everything asked for is present. Jaccard, which counts the union
# in the denominator, prefers Shortbread: Sachertorte's seven extra ingredients
# are seven things the two recipes do not share, and Jaccard charges for them.
#
# Neither is wrong. They answer different questions -- "is what I asked for
# there?" against "how much of everything involved is shared?" -- and the day's
# job is to notice that you have to pick one.

RECIPE_QUERY = frozenset({"flour", "butter", "sugar", "egg"})
RECIPES: dict[str, frozenset[str]] = {
    "Sachertorte": frozenset({
        "flour", "butter", "sugar", "egg",
        "cocoa", "apricot jam", "chocolate", "vanilla",
        "salt", "milk", "almond",
    }),
    "Shortbread": frozenset({"flour", "butter", "cornflour"}),
}

# ---------------------------------------------------------------------------
# 5. Mahalanobis: Euclidean after accounting for how the data actually varies.
# ---------------------------------------------------------------------------
#
# Eight readings of two sensors that move together almost perfectly. The mean
# is exactly (0, 0) and the population covariance comes out exactly
#
#     [[7.5, 7.0],
#      [7.0, 7.5]]
#
# whose determinant is exactly 7.25. The data lies along the line y = x: that
# is the grain of it, and Day 106's eigenvectors of this matrix are what name
# that direction.
#
# The two probe points are the same Euclidean distance from the mean -- both
# sqrt(18) = 4.2426... -- and nothing about ordinary distance can tell them
# apart. Mahalanobis can: ALONG the grain is cheap, ACROSS it is expensive.

SENSOR_READINGS: tuple[tuple[float, float], ...] = (
    (-4.0, -3.0),
    (-3.0, -4.0),
    (-2.0, -1.0),
    (-1.0, -2.0),
    (1.0, 2.0),
    (2.0, 1.0),
    (3.0, 4.0),
    (4.0, 3.0),
)

# Along the grain of the data: both sensors high together, which is what this
# pair of sensors does all day.
PROBE_ALONG = (3.0, 3.0)
# Across the grain: one sensor high while the other is low, which never happens
# in the eight readings above. Same Euclidean distance. Not the same event.
PROBE_ACROSS = (3.0, -3.0)

# ---------------------------------------------------------------------------
# 6. The scaling demonstration: the thing that silently decides your answer.
# ---------------------------------------------------------------------------
#
# A bearing catalogue with two features recorded in the units the supplier
# happened to use: bore diameter in METRES and mass in GRAMS. The numbers in
# one column are around 0.02 and in the other around 350, so squared
# differences in the second column are roughly ten million times larger. The
# bore column does not lose the argument. It never enters it.

BEARING_FEATURES = ("bore diameter (m)", "mass (g)")

BEARING_QUERY = (0.020, 300.0)

BEARINGS: dict[str, tuple[float, float]] = {
    # Bore matches the query EXACTLY. 40 g heavier.
    "P": (0.020, 340.0),
    # Bore is 12 mm too big -- 60 per cent out, and unusable. Mass is 2 g off.
    "R": (0.032, 302.0),
    "S": (0.008, 250.0),
    "T": (0.026, 410.0),
    "U": (0.014, 275.0),
    "V": (0.038, 500.0),
}

# ---------------------------------------------------------------------------
# 7. The counter-example that shows cosine distance is not a metric.
# ---------------------------------------------------------------------------
#
# Day 103 proved this. It is restated here rather than re-derived, because a
# concrete triple is worth more than the proof once you have seen the proof.
#
#   cosine_distance(EAST, DIAGONAL) + cosine_distance(DIAGONAL, NORTH)
#     = 0.2929 + 0.2929 = 0.5858
#   cosine_distance(EAST, NORTH)
#     = 1.0
#
# The direct route is longer than going via a third point, which no metric may
# ever allow.

EAST = (1.0, 0.0)
DIAGONAL = (1.0, 1.0)
NORTH = (0.0, 1.0)

# The triple used for the POSITIVE side of the same check: L1, L2 and
# L-infinity must all satisfy the triangle inequality on every triple, and
# these three vectors are checked exhaustively in all six orderings.
TRIANGLE_TRIPLE = ((1.0, 7.0, 2.0), (4.0, 1.0, 9.0), (-2.0, 3.0, 3.0))

# The single vector every norm axiom is checked on, and the scalar it is
# multiplied by for absolute homogeneity.
AXIOM_VECTOR = (3.0, -4.0, 12.0)
AXIOM_SCALAR = -2.5

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

Both `examples/` and `starter/` contain modules called `measures` and
`catalogue`, 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 `measures` was seen first and then reuse it for the
other suite -- so these starter tests would silently pass against the reference
solution instead of skipping. That is a wrong answer with a green tick on it,
which is the worst kind.

So: put this directory first on the import path, and drop any already-imported
`measures`, `catalogue` or `answers` 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 ("measures", "catalogue", "answers"):
    module = sys.modules.get(name)
    origin = getattr(module, "__file__", "") or ""
    if module is not None and not origin.startswith(HERE):
        del sys.modules[name]

"""Exercise 1 -- every measure in this lab, built from arithmetic alone.

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

Check yourself as you go, from the LAB DIRECTORY (the one above this file):

    .venv/bin/pytest starter -q

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

**Use only the standard library here.** `math`, `abs`, `sum`, `max`, `sorted`
and set operations are everything you need. NumPy is not forbidden by a lint
rule -- it is forbidden by the point of the exercise. The tests check your work
against `numpy.linalg.norm(v, ord=p)`, `numpy.cov` and `numpy.linalg.inv`, and
that check means nothing if your answer was NumPy's answer all along.

Two conventions, fixed and used everywhere:

1. A vector is a plain sequence of floats. Two vectors of different lengths may
   not be compared, and `_paired` below raises rather than zipping the shorter
   one and quietly answering a different question.

2. A DISTANCE gets smaller as things get more alike; a SIMILARITY gets larger.
   Every function name says which it is. Nothing in this module guesses.
"""

from __future__ import annotations

import math
from collections.abc import Callable, Hashable, Iterable, Sequence

Vector = Sequence[float]
Matrix = list[list[float]]

# Every float comparison in this lab is made against this tolerance. Written
# for you, and it is not decoration: the Mahalanobis result in exercise 1.16
# comes out as exactly 6.0 through one correct route and 5.999999999999999
# through another. `== 6.0` would pass for one and fail for the other.
TOL = 1e-12


class DimensionMismatch(ValueError):
    """Raised when two vectors of different lengths are compared.

    Written for you. Subclasses ValueError so that an existing
    `except ValueError` catches it, matching how NumPy reports the same
    mistake.
    """


def _paired(u: Vector, v: Vector) -> list[tuple[float, float]]:
    """Pair two vectors elementwise, refusing to compare different lengths.

    Written for you. Use it -- `zip` alone will silently truncate.
    """
    if len(u) != len(v):
        raise DimensionMismatch(
            f"vectors have different lengths: {len(u)} and {len(v)}"
        )
    return list(zip(u, v))


# -- Exercise 1.1 to 1.4: the norms -------------------------------------------


def l1_norm(v: Vector) -> float:
    """The L1 norm of `v`: the sum of the absolute values.

        l1_norm((3.0, -4.0, 12.0))  ->  19.0

    One line with `sum` and a generator expression.
    """
    raise NotImplementedError("exercise 1.1: l1_norm")


def l2_norm(v: Vector) -> float:
    """The L2 norm of `v`: the square root of the sum of squares.

        l2_norm((3.0, -4.0, 12.0))  ->  13.0     exactly, since 169 = 13 * 13

    Use `math.sqrt`. Do not use `** 0.5`; they agree here but `math.sqrt` says
    what you meant.
    """
    raise NotImplementedError("exercise 1.2: l2_norm")


def linf_norm(v: Vector) -> float:
    """The L-infinity norm of `v`: the largest single absolute component.

        linf_norm((3.0, -4.0, 12.0))  ->  12.0
        linf_norm(())                 ->  0.0

    Note the empty case. `max` on an empty sequence raises, so pass
    `default=0.0`.
    """
    raise NotImplementedError("exercise 1.3: linf_norm")


def p_norm(v: Vector, p: float) -> float:
    """The general p-norm: (sum of |x| ** p) ** (1 / p).

        p_norm((3, 4), 1)        ->  7.0
        p_norm((3, 4), 2)        ->  5.0
        p_norm((3, 4), math.inf) ->  4.0

    Two cases the tests check and one trap:

    * `p < 1` is NOT a norm -- the triangle inequality fails below 1 -- so
      raise `ValueError` rather than returning a plausible number.
    * `p = math.inf` must be handled as the LIMIT, returning the largest
      absolute component. Computing it as arithmetic overflows, because
      `x ** math.inf` is `inf` for any x above 1.
    * `math.isinf(p)` is the readable test for the second case.
    """
    raise NotImplementedError("exercise 1.4: p_norm")


# -- Exercise 1.5 to 1.7: the distances ---------------------------------------


def l1_distance(u: Vector, v: Vector) -> float:
    """Manhattan distance: total disagreement, summed across the features.

        l1_distance((4, 3, 2, 1), (4, 3, 2, 6))  ->  5.0

    Every unit of difference costs the same wherever it happens. Use `_paired`.
    """
    raise NotImplementedError("exercise 1.5: l1_distance")


def l2_distance(u: Vector, v: Vector) -> float:
    """Euclidean distance: straight-line separation.

        l2_distance((0, 0), (6, 8))  ->  10.0

    Squaring makes one large disagreement cost far more than several small
    ones that add up to the same total.
    """
    raise NotImplementedError("exercise 1.6: l2_distance")


def linf_distance(u: Vector, v: Vector) -> float:
    """Chebyshev distance: the single worst feature decides, alone.

        linf_distance((0, 0), (6, 8))  ->  8.0

    Every other feature is ignored entirely. Remember `default=0.0`.
    """
    raise NotImplementedError("exercise 1.7: linf_distance")


def minkowski_distance(u: Vector, v: Vector, p: float) -> float:
    """The p-norm of the difference. Written for you, once p_norm exists."""
    return p_norm([a - b for a, b in _paired(u, v)], p)


# -- Exercise 1.8: the angle --------------------------------------------------


def dot(u: Vector, v: Vector) -> float:
    """Day 103's dot product. Written for you."""
    return sum(a * b for a, b in _paired(u, v))


def cosine_similarity(u: Vector, v: Vector) -> float:
    """The cosine of the angle between two vectors: 1 identical, 0 orthogonal.

        cosine_similarity((4, 3, 2, 1), (12, 9, 6, 3))  ->  1.0

    Day 103 derived this: the dot product divided by both lengths. Length is
    divided out, which is the whole point -- the second vector above is the
    first at three times the size and scores exactly 1.

    The zero vector has no direction, so cosine is undefined for it. RAISE
    `ValueError` rather than returning 0.0; a silent wrong answer here
    propagates into a ranking and is very hard to find later. Use `TOL` to
    test for a zero length rather than `== 0`.
    """
    raise NotImplementedError("exercise 1.8: cosine_similarity")


def cosine_distance(u: Vector, v: Vector) -> float:
    """1 minus the cosine similarity. Written for you.

    Widely used, useful, and NOT a metric: exercise 2 asks you to say which
    axiom it breaks.
    """
    return 1.0 - cosine_similarity(u, v)


# -- Exercise 1.9 to 1.11: categorical and set data ---------------------------


def hamming_distance(a: Sequence, b: Sequence) -> int:
    """How many positions differ. The right answer for categorical features.

        hamming_distance(("steel", "zinc", "M8"),
                         ("brass", "zinc", "M8"))   ->  1

    Nothing is subtracted, so the values need not be numbers. Return an `int`,
    not a float, and use `_paired` so that a length mismatch raises.
    """
    raise NotImplementedError("exercise 1.9: hamming_distance")


def normalised_hamming(a: Sequence, b: Sequence) -> float:
    """Hamming as a fraction of the fields. Written for you."""
    if not a:
        raise ValueError("normalised Hamming needs at least one field")
    return hamming_distance(a, b) / len(a)


def jaccard_similarity(a: Iterable[Hashable], b: Iterable[Hashable]) -> float:
    """|intersection| / |union| for two sets: 1 identical, 0 disjoint.

        jaccard_similarity({1, 2, 3, 4}, {1, 2, 3, 4, 5})  ->  0.8
        jaccard_similarity({1, 2}, {3, 4})                 ->  0.0

    Two conventions the tests check:

    * The arguments may be any iterables. Call `set()` on both first.
    * Two EMPTY sets are defined here as identical, similarity 1.0. That is a
      convention rather than a derivation -- 0/0 has no answer -- and stating
      it is better than letting a ZeroDivisionError escape at 3 a.m.
    """
    raise NotImplementedError("exercise 1.10: jaccard_similarity")


def jaccard_distance(a: Iterable[Hashable], b: Iterable[Hashable]) -> float:
    """1 minus Jaccard similarity. Written for you.

    Unlike cosine distance, this one IS a metric -- the reference tests check
    that on all 4096 triples of subsets of a four-element set.
    """
    return 1.0 - jaccard_similarity(a, b)


def vocabulary(*collections: Iterable[Hashable]) -> list[Hashable]:
    """The SORTED union of several collections. Written for you.

    Sorted rather than in encounter order, so the binary vectors below are the
    same on every run. A set has no order, and a vector built from one without
    sorting is a different vector each time the interpreter starts.
    """
    seen: set[Hashable] = set()
    for collection in collections:
        seen |= set(collection)
    return sorted(seen)


def to_binary_vector(items: Iterable[Hashable],
                     axes: Sequence[Hashable]) -> list[float]:
    """Turn a set into a 1/0 vector over a fixed list of axes.

        to_binary_vector({"a", "c"}, ["a", "b", "c"])  ->  [1.0, 0.0, 1.0]

    This is how a set gets handed to a measure that expects numbers, and it is
    what makes the Jaccard-against-cosine comparison possible on identical
    data. Return floats, not booleans.
    """
    raise NotImplementedError("exercise 1.11: to_binary_vector")


# -- Written for you: small matrix arithmetic ---------------------------------
#
# You have built matrix multiplication already, on Day 101, and Gauss-Jordan
# elimination is a day of its own. These are given so that exercise 1.16 is
# about Mahalanobis distance and not about linear solvers.


def transpose(m: Matrix) -> Matrix:
    """Rows become columns."""
    return [list(col) for col in zip(*m)]


def matmul(a: Matrix, b: Matrix) -> Matrix:
    """Day 101's matrix product."""
    if len(a[0]) != len(b):
        raise DimensionMismatch(
            f"cannot multiply {len(a)}x{len(a[0])} by {len(b)}x{len(b[0])}"
        )
    bt = transpose(b)
    return [[sum(x * y for x, y in zip(row, col)) for col in bt] for row in a]


def mat_vec(m: Matrix, v: Vector) -> list[float]:
    """Matrix times column vector."""
    if len(m[0]) != len(v):
        raise DimensionMismatch(
            f"cannot apply {len(m)}x{len(m[0])} matrix to a vector of {len(v)}"
        )
    return [sum(x * y for x, y in zip(row, v)) for row in m]


def inverse(m: Matrix) -> Matrix:
    """Invert a square matrix by Gauss-Jordan elimination with partial pivoting.

    Raises ValueError when the matrix is singular, which is what a covariance
    matrix with two identical features gives you.
    """
    n = len(m)
    if any(len(row) != n for row in m):
        raise DimensionMismatch("only a square matrix can be inverted")
    aug = [list(map(float, row)) + [1.0 if i == j else 0.0 for j in range(n)]
           for i, row in enumerate(m)]
    for col in range(n):
        pivot = max(range(col, n), key=lambda r: abs(aug[r][col]))
        if abs(aug[pivot][col]) <= TOL:
            raise ValueError("matrix is singular: it has no inverse")
        aug[col], aug[pivot] = aug[pivot], aug[col]
        scale = aug[col][col]
        aug[col] = [x / scale for x in aug[col]]
        for row in range(n):
            if row == col:
                continue
            factor = aug[row][col]
            if factor == 0.0:
                continue
            aug[row] = [x - factor * y for x, y in zip(aug[row], aug[col])]
    return [row[n:] for row in aug]


# -- Exercise 1.12 to 1.16: statistics, scaling and Mahalanobis --------------


def column_means(rows: Sequence[Vector]) -> list[float]:
    """The mean of each column of a table of equal-length rows.

        column_means([(1.0, 10.0), (3.0, 20.0)])  ->  [2.0, 15.0]

    Raise `ValueError` on an empty table rather than returning `[]`.
    """
    raise NotImplementedError("exercise 1.12: column_means")


def column_stds(rows: Sequence[Vector]) -> list[float]:
    """The POPULATION standard deviation of each column: divided by n.

        column_stds([(1.0,), (3.0,)])  ->  [1.0]

    Which divisor to use is a real decision and this lab makes it explicitly.
    `n` describes the table you have; `n - 1` estimates a wider population you
    are sampling from. Scikit-learn's StandardScaler divides by `n`, and
    `numpy.std` does too unless you pass `ddof=1`, so `n` is what this lab
    uses and what the tests check against.
    """
    raise NotImplementedError("exercise 1.13: column_stds")


def standardise(
    rows: Sequence[Vector],
    means: Sequence[float] | None = None,
    stds: Sequence[float] | None = None,
) -> list[list[float]]:
    """Subtract the column mean, divide by the column standard deviation.

        standardise([(1.0,), (3.0,)])  ->  [[-1.0], [1.0]]

    Also called the z-score. Three requirements the tests check:

    * When `means` or `stds` is given, USE IT instead of recomputing. This is
      how a query gets standardised with the same numbers as the catalogue it
      is being compared against. Standardising a single query against itself
      gives a row of zeros, which is a real bug with a long history.
    * A column whose standard deviation is 0 (within `TOL`) must come out as
      0.0 rather than dividing by zero. It carries no information to scale.
    * Return a list of lists of floats, one per input row, same order.
    """
    raise NotImplementedError("exercise 1.14: standardise")


def covariance_matrix(rows: Sequence[Vector]) -> Matrix:
    """The population covariance matrix of a table of rows, divided by n.

    Entry (i, j) is the average product of column i's and column j's
    deviations from their own means:

        cov[i][j] = sum over rows of (row[i] - mean[i]) * (row[j] - mean[j]) / n

    The diagonal is each column's variance. On the eight sensor readings in
    `catalogue.py` the answer is exactly [[7.5, 7.0], [7.0, 7.5]], which you
    can check by hand -- and the tests do.

    Day 106's eigenvectors of this matrix are the directions the data actually
    spreads along, which is exactly what the next function measures in.
    """
    raise NotImplementedError("exercise 1.15: covariance_matrix")


def mahalanobis_distance(u: Vector, v: Vector, cov_inverse: Matrix) -> float:
    """Euclidean distance after accounting for how the data actually varies.

    One line of arithmetic. Take the difference `z = u - v`, and instead of
    dotting it with itself, dot it with itself THROUGH the inverse covariance:

        d = sqrt( z . (cov_inverse . z) )

    You have `dot` and `mat_vec` already.

    Check yourself two ways:

    * Passing the IDENTITY matrix as `cov_inverse` must give back ordinary
      Euclidean distance, exactly. That is the cleanest statement of what the
      covariance is doing.
    * On the sensor readings, (3, 3) and (3, -3) are the same Euclidean
      distance from the mean and 1.114172 against 6.0 in Mahalanobis.

    One floating-point guard the tests check. A covariance matrix is positive
    semi-definite, so the value under the square root is non-negative in exact
    arithmetic -- but a result that should be 0 can come out as -1e-17, and
    `math.sqrt` raises on it. Clamp a tiny negative to 0.0, and raise
    `ValueError` for one that is genuinely negative (worse than -TOL), because
    that means the matrix you were handed is not an inverse covariance.
    """
    raise NotImplementedError("exercise 1.16: mahalanobis_distance")


# -- Exercise 1.17: the function the whole day is about ----------------------


def rank(
    query,
    candidates: dict[str, object],
    measure: Callable,
    higher_is_better: bool = False,
) -> list[tuple[str, float]]:
    """Score every candidate against the query and sort best first.

        rank((4, 3, 2, 1), {"Aisle": (4, 3, 2, 6)}, l1_distance)
            ->  [("Aisle", 5.0)]

    This is the function the day is really about. Swapping Manhattan for
    cosine must be ONE argument here.

    Three requirements the tests check:

    * Return a list of `(name, score)` pairs, score as a `float`.
    * `higher_is_better=False` sorts ascending (a distance); `True` sorts
      descending (a similarity). Getting this backwards returns the WORST
      match with complete confidence and no error message.
    * Ties must break by the candidate's name, so two runs never disagree for
      a reason that has nothing to do with the data. A one-key sort on the
      score alone is not enough; sort on a tuple.
    """
    raise NotImplementedError("exercise 1.17: rank")


def winner(query, candidates: dict[str, object], measure: Callable,
           higher_is_better: bool = False) -> str:
    """The name at the top of `rank`. Written for you."""
    return rank(query, candidates, measure, higher_is_better)[0][0]

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

    .venv/bin/pytest starter -q

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

Float comparisons use the tolerance TOL stated in measures.py, except the
predictions in answers.py, which use 1e-9 so you need not type more decimals
than the question asks for. Counts, names and rankings are compared exactly,
because they are exact.
"""

from __future__ import annotations

import itertools
import math

import numpy as np
import pytest

import answers
import catalogue
import measures

TOL = measures.TOL
PREDICTION_TOL = 1e-9

Q = catalogue.QUERY
ARTICLES = catalogue.ARTICLES


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


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


def close(a, b, tol=PREDICTION_TOL):
    """Elementwise closeness for numbers, vectors and matrices."""
    if isinstance(a, (list, tuple)) and isinstance(b, (list, tuple)):
        return len(a) == len(b) and all(close(x, y, tol) for x, y in zip(a, b))
    return abs(a - b) <= tol


def ranked(query, candidates, measure, higher_is_better=False):
    """`rank`, routed through `written` so an unwritten `rank` skips."""
    return written(measures.rank, query, candidates, measure, higher_is_better)


def top(query, candidates, measure, higher_is_better=False):
    return ranked(query, candidates, measure, higher_is_better)[0][0]


# -- 0. Always passes: the data is what the lab says it is --------------------


def test_0_00_the_catalogue_is_intact():
    """One test that passes before you write anything, so a green run means
    the suite itself is working rather than that nothing was collected."""
    assert catalogue.QUERY == (4, 3, 2, 1)
    assert catalogue.ARTICLES["Cartogram"] == tuple(3 * c for c in Q)
    assert len(catalogue.SENSOR_READINGS) == 8
    assert len(catalogue.BEARINGS) == 6
    assert len(catalogue.RECIPES["Sachertorte"]) == 11


# -- 1. The norms -------------------------------------------------------------


def test_1_01_l1_norm():
    assert written(measures.l1_norm, (3.0, -4.0, 12.0)) == 19.0
    assert written(measures.l1_norm, ()) == 0


def test_1_02_l2_norm():
    assert close(written(measures.l2_norm, (3.0, -4.0, 12.0)), 13.0, TOL)
    assert close(written(measures.l2_norm, (0.0, 0.0)), 0.0, TOL)


def test_1_03_linf_norm():
    assert written(measures.linf_norm, (3.0, -4.0, 12.0)) == 12.0
    assert written(measures.linf_norm, ()) == 0.0


def test_1_04_p_norm_reproduces_the_three_named_norms():
    v = (3.0, 4.0)
    assert close(written(measures.p_norm, v, 1), 7.0, TOL)
    assert close(written(measures.p_norm, v, 2), 5.0, TOL)
    assert close(written(measures.p_norm, v, math.inf), 4.0, TOL)


def test_1_05_p_norm_matches_numpy_for_fractional_p():
    v = (3.0, 4.0)
    for p in (1.5, 3, 8):
        mine = written(measures.p_norm, v, p)
        assert close(mine, float(np.linalg.norm(np.asarray(v), ord=p)), TOL)


def test_1_06_p_norm_refuses_p_below_one():
    try:
        measures.p_norm((3.0, 4.0), 0.5)
    except NotImplementedError as exc:
        pytest.skip(f"not written yet: {exc}")
    except ValueError:
        return
    pytest.fail("p_norm(v, 0.5) should raise ValueError: below p = 1 it is "
                "not a norm, because the triangle inequality fails")


def test_1_07_p_norm_is_non_increasing_in_p():
    v = (3.0, 4.0)
    values = [written(measures.p_norm, v, p) for p in (1, 1.5, 2, 3, 8, 64)]
    for earlier, later in zip(values, values[1:]):
        assert later <= earlier + TOL


# -- 2. The distances ---------------------------------------------------------


def test_2_01_l1_distance():
    assert close(written(measures.l1_distance, Q, ARTICLES["Aisle"]), 5.0, TOL)
    assert close(written(measures.l1_distance, Q, ARTICLES["Beacon"]), 6.0, TOL)


def test_2_02_l2_distance():
    assert close(written(measures.l2_distance, (0.0, 0.0), (6.0, 8.0)),
                 10.0, TOL)
    assert close(written(measures.l2_distance, Q, ARTICLES["Beacon"]),
                 math.sqrt(12.0), TOL)


def test_2_03_linf_distance():
    assert close(written(measures.linf_distance, (0.0, 0.0), (6.0, 8.0)),
                 8.0, TOL)
    assert close(written(measures.linf_distance, Q, ARTICLES["Aisle"]),
                 5.0, TOL)


def test_2_04_all_three_match_numpy_on_the_articles():
    for vec in ARTICLES.values():
        d = np.asarray(Q, dtype=float) - np.asarray(vec, dtype=float)
        assert close(written(measures.l1_distance, Q, vec),
                     float(np.linalg.norm(d, ord=1)), TOL)
        assert close(written(measures.l2_distance, Q, vec),
                     float(np.linalg.norm(d, ord=2)), TOL)
        assert close(written(measures.linf_distance, Q, vec),
                     float(np.linalg.norm(d, ord=np.inf)), TOL)


def test_2_05_comparing_different_lengths_raises():
    try:
        measures.l2_distance((1.0, 2.0, 3.0), (1.0, 2.0))
    except NotImplementedError as exc:
        pytest.skip(f"not written yet: {exc}")
    except measures.DimensionMismatch:
        return
    pytest.fail("comparing a length-3 vector with a length-2 vector must "
                "raise DimensionMismatch: use _paired rather than zip")


def test_2_06_the_ordering_linf_le_l2_le_l1_holds():
    for vec in ARTICLES.values():
        a = written(measures.linf_distance, Q, vec)
        b = written(measures.l2_distance, Q, vec)
        c = written(measures.l1_distance, Q, vec)
        assert a <= b + TOL <= c + 2 * TOL


# -- 3. Cosine ---------------------------------------------------------------


def test_3_01_cosine_of_a_scaled_copy_is_one():
    assert close(written(measures.cosine_similarity, Q, ARTICLES["Cartogram"]),
                 1.0, TOL)


def test_3_02_cosine_matches_numpy():
    for vec in ARTICLES.values():
        a = np.asarray(Q, dtype=float)
        b = np.asarray(vec, dtype=float)
        theirs = float(a @ b / (np.linalg.norm(a) * np.linalg.norm(b)))
        assert close(written(measures.cosine_similarity, Q, vec), theirs, TOL)


def test_3_03_cosine_of_orthogonal_vectors_is_zero():
    assert close(written(measures.cosine_similarity, (1.0, 0.0), (0.0, 1.0)),
                 0.0, TOL)


def test_3_04_cosine_refuses_the_zero_vector():
    try:
        measures.cosine_similarity((0.0, 0.0), (1.0, 1.0))
    except NotImplementedError as exc:
        pytest.skip(f"not written yet: {exc}")
    except ValueError:
        return
    pytest.fail("cosine_similarity must raise ValueError for a zero vector: "
                "it has no direction, and returning 0.0 hides the problem")


def test_3_05_cosine_distance_breaks_the_triangle_inequality():
    east, diagonal, north = (catalogue.EAST, catalogue.DIAGONAL,
                             catalogue.NORTH)
    detour = (written(measures.cosine_distance, east, diagonal)
              + written(measures.cosine_distance, diagonal, north))
    direct = written(measures.cosine_distance, east, north)
    assert detour < direct - TOL
    assert close(direct, 1.0, TOL)


# -- 4. Categorical and set data ---------------------------------------------


def test_4_01_hamming_on_the_parts_register():
    ref = catalogue.REFERENCE_RECORD
    got = {n: written(measures.hamming_distance, ref, r)
           for n, r in catalogue.CANDIDATE_RECORDS.items()}
    assert got == {"part-71": 1, "part-72": 3, "part-73": 6}


def test_4_02_hamming_returns_an_int_and_needs_no_arithmetic():
    value = written(measures.hamming_distance, ("red", "blue"),
                    ("green", "blue"))
    assert value == 1
    assert isinstance(value, int)


def test_4_03_hamming_on_the_bit_flags():
    assert written(measures.hamming_distance, catalogue.FLAGS_A,
                   catalogue.FLAGS_B) == 2


def test_4_04_jaccard_basic_cases():
    assert close(written(measures.jaccard_similarity, {1, 2, 3, 4},
                         {1, 2, 3, 4, 5}), 0.8, TOL)
    assert close(written(measures.jaccard_similarity, {1, 2}, {3, 4}),
                 0.0, TOL)
    assert close(written(measures.jaccard_similarity, set(), set()), 1.0, TOL)


def test_4_05_jaccard_on_the_recipes():
    q = catalogue.RECIPE_QUERY
    assert close(written(measures.jaccard_similarity, q,
                         catalogue.RECIPES["Sachertorte"]), 4 / 11, TOL)
    assert close(written(measures.jaccard_similarity, q,
                         catalogue.RECIPES["Shortbread"]), 2 / 5, TOL)


def test_4_06_to_binary_vector():
    assert written(measures.to_binary_vector, {"a", "c"},
                   ["a", "b", "c"]) == [1.0, 0.0, 1.0]


def test_4_07_jaccard_and_cosine_disagree_on_the_recipes():
    q = catalogue.RECIPE_QUERY
    axes = measures.vocabulary(q, *catalogue.RECIPES.values())
    qv = written(measures.to_binary_vector, q, axes)
    jac = {n: written(measures.jaccard_similarity, q, s)
           for n, s in catalogue.RECIPES.items()}
    cos = {n: written(measures.cosine_similarity, qv,
                      written(measures.to_binary_vector, s, axes))
           for n, s in catalogue.RECIPES.items()}
    assert max(jac, key=jac.get) == "Shortbread"
    assert max(cos, key=cos.get) == "Sachertorte"


def test_4_08_jaccard_distance_satisfies_the_triangle_inequality():
    universe = ("a", "b", "c")
    subsets = [frozenset(c) for r in range(len(universe) + 1)
               for c in itertools.combinations(universe, r)]
    for a, b, c in itertools.product(subsets, repeat=3):
        assert (written(measures.jaccard_distance, a, b)
                + written(measures.jaccard_distance, b, c)
                >= written(measures.jaccard_distance, a, c) - TOL)


# -- 5. Statistics and scaling -----------------------------------------------


def test_5_01_column_means():
    assert close(written(measures.column_means,
                         [(1.0, 10.0), (3.0, 20.0)]), [2.0, 15.0], TOL)


def test_5_02_column_stds_use_the_population_divisor():
    rows = list(catalogue.BEARINGS.values())
    mine = written(measures.column_stds, rows)
    theirs = np.asarray(rows, dtype=float).std(axis=0)
    assert close(list(mine), list(theirs), TOL)


def test_5_03_column_stds_are_not_the_sample_divisor():
    rows = list(catalogue.BEARINGS.values())
    mine = written(measures.column_stds, rows)
    sample = np.asarray(rows, dtype=float).std(axis=0, ddof=1)
    assert not close(list(mine), list(sample), 1e-6)


def test_5_04_standardise_gives_mean_zero_and_sd_one():
    rows = list(catalogue.BEARINGS.values())
    z = written(measures.standardise, rows)
    assert close(written(measures.column_means, z), [0.0, 0.0], 1e-12)
    assert close(written(measures.column_stds, z), [1.0, 1.0], 1e-12)


def test_5_05_standardise_uses_supplied_means_and_stds():
    rows = [(1.0,), (3.0,)]
    z = written(measures.standardise, rows, [0.0], [1.0])
    assert close(z, [[1.0], [3.0]], TOL)


def test_5_06_standardise_leaves_a_constant_column_at_zero():
    z = written(measures.standardise, [(1.0, 5.0), (2.0, 5.0), (3.0, 5.0)])
    assert [row[1] for row in z] == [0.0, 0.0, 0.0]


def test_5_07_standardising_changes_the_bearing_winner():
    rows = list(catalogue.BEARINGS.values())
    means = written(measures.column_means, rows)
    stds = written(measures.column_stds, rows)
    raw = top(catalogue.BEARING_QUERY, catalogue.BEARINGS,
              measures.l2_distance)
    q = written(measures.standardise, [catalogue.BEARING_QUERY],
                means, stds)[0]
    scaled = {n: written(measures.standardise, [v], means, stds)[0]
              for n, v in catalogue.BEARINGS.items()}
    assert raw == "R"
    assert top(q, scaled, measures.l2_distance) == "P"


# -- 6. Covariance and Mahalanobis -------------------------------------------


def test_6_01_covariance_of_the_sensor_readings():
    cov = written(measures.covariance_matrix, catalogue.SENSOR_READINGS)
    assert close(cov, [[7.5, 7.0], [7.0, 7.5]], TOL)


def test_6_02_covariance_matches_numpy_with_bias_true():
    mine = written(measures.covariance_matrix, catalogue.SENSOR_READINGS)
    theirs = np.cov(np.asarray(catalogue.SENSOR_READINGS, dtype=float),
                    rowvar=False, bias=True)
    assert close(mine, theirs.tolist(), TOL)


def test_6_03_mahalanobis_with_the_identity_is_euclidean():
    identity = [[1.0, 0.0], [0.0, 1.0]]
    for probe in ((3.0, 3.0), (3.0, -3.0), (-2.5, 4.75)):
        assert close(
            written(measures.mahalanobis_distance, probe, (0.0, 0.0),
                    identity),
            written(measures.l2_distance, probe, (0.0, 0.0)), TOL)


def test_6_04_mahalanobis_separates_the_two_probes():
    mean = written(measures.column_means, catalogue.SENSOR_READINGS)
    cov = written(measures.covariance_matrix, catalogue.SENSOR_READINGS)
    inv = measures.inverse(cov)
    along = written(measures.mahalanobis_distance, catalogue.PROBE_ALONG,
                    mean, inv)
    across = written(measures.mahalanobis_distance, catalogue.PROBE_ACROSS,
                     mean, inv)
    assert close(across, 6.0, TOL)
    assert close(along, math.sqrt(9.0 / 7.25), TOL)


def test_6_05_euclidean_cannot_separate_them():
    mean = written(measures.column_means, catalogue.SENSOR_READINGS)
    assert close(written(measures.l2_distance, catalogue.PROBE_ALONG, mean),
                 written(measures.l2_distance, catalogue.PROBE_ACROSS, mean),
                 TOL)


def test_6_06_mahalanobis_clamps_a_tiny_negative_rather_than_raising():
    inv = measures.inverse(written(measures.covariance_matrix,
                                   catalogue.SENSOR_READINGS))
    assert close(written(measures.mahalanobis_distance, (1.0, 1.0),
                         (1.0, 1.0), inv), 0.0, TOL)


# -- 7. The ranking function -------------------------------------------------


def test_7_01_rank_returns_name_score_pairs_for_every_candidate():
    order = ranked(Q, ARTICLES, measures.l1_distance)
    assert sorted(n for n, _ in order) == sorted(ARTICLES)
    assert all(isinstance(score, float) for _, score in order)


def test_7_02_rank_ascends_for_a_distance():
    scores = [s for _, s in ranked(Q, ARTICLES, measures.l2_distance)]
    assert scores == sorted(scores)


def test_7_03_rank_descends_for_a_similarity():
    scores = [s for _, s in ranked(Q, ARTICLES, measures.cosine_similarity,
                                   higher_is_better=True)]
    assert scores == sorted(scores, reverse=True)


def test_7_04_ties_break_by_name():
    candidates = {"zulu": (1.0, 1.0), "alpha": (1.0, 1.0), "mike": (1.0, 1.0)}
    order = [n for n, _ in ranked((0.0, 0.0), candidates,
                                  measures.l2_distance)]
    assert order == ["alpha", "mike", "zulu"]


def test_7_05_three_measures_name_three_different_winners():
    picks = {
        "l1": top(Q, ARTICLES, measures.l1_distance),
        "l2": top(Q, ARTICLES, measures.l2_distance),
        "cos": top(Q, ARTICLES, measures.cosine_similarity, True),
    }
    assert picks == {"l1": "Aisle", "l2": "Beacon", "cos": "Cartogram"}


def test_7_06_the_warehouse_displacement():
    a, b = catalogue.FLOOR_FROM, catalogue.FLOOR_TO
    assert close(written(measures.l1_distance, a, b), 14.0, TOL)
    assert close(written(measures.l2_distance, a, b), 10.0, TOL)
    assert close(written(measures.linf_distance, a, b), 8.0, TOL)


def test_7_07_chebyshev_accepts_the_part_the_others_would_reject():
    nominal = catalogue.NOMINAL_PART
    limit = catalogue.PART_TOLERANCE_MM
    a = catalogue.MEASURED_PARTS["batch-A"]
    b = catalogue.MEASURED_PARTS["batch-B"]
    assert written(measures.linf_distance, a, nominal) <= limit + TOL
    assert written(measures.linf_distance, b, nominal) > limit
    assert (written(measures.l1_distance, b, nominal)
            < written(measures.l1_distance, a, nominal))


# -- 8. Your predictions ------------------------------------------------------


def test_8_01_l1_query_to_aisle():
    assert close(predicted("L1_QUERY_TO_AISLE"), 5.0)


def test_8_02_l1_query_to_beacon():
    assert close(predicted("L1_QUERY_TO_BEACON"), 6.0)


def test_8_03_l2_query_to_aisle():
    assert close(predicted("L2_QUERY_TO_AISLE"), 5.0)


def test_8_04_linf_query_to_beacon():
    assert close(predicted("LINF_QUERY_TO_BEACON"), 2.0)


def test_8_05_cosine_query_to_cartogram():
    assert close(predicted("COSINE_QUERY_TO_CARTOGRAM"), 1.0, 5e-5)


def test_8_06_l1_winner():
    assert predicted("L1_WINNER") == "Aisle"


def test_8_07_l2_winner():
    assert predicted("L2_WINNER") == "Beacon"


def test_8_08_cosine_winner():
    assert predicted("COSINE_WINNER") == "Cartogram"


def test_8_09_p_norm_at_1():
    assert close(predicted("P_NORM_3_4_AT_1"), 7.0)


def test_8_10_p_norm_at_2():
    assert close(predicted("P_NORM_3_4_AT_2"), 5.0)


def test_8_11_p_norm_at_infinity():
    assert close(predicted("P_NORM_3_4_AT_INF"), 4.0)


def test_8_12_p_norm_falls_as_p_rises():
    assert predicted("P_NORM_AS_P_RISES") == "fall"


def test_8_13_squared_euclidean_breaks_absolute_homogeneity():
    assert predicted("AXIOM_SQUARED_EUCLIDEAN_BREAKS") == "absolute homogeneity"


def test_8_14_cosine_distance_east_north():
    assert close(predicted("COSINE_DISTANCE_EAST_NORTH"), 1.0)


def test_8_15_cosine_triangle_does_not_hold():
    assert predicted("COSINE_TRIANGLE_HOLDS") is False


def test_8_16_jaccard_distance_is_a_metric():
    assert predicted("JACCARD_DISTANCE_IS_A_METRIC") is True


def test_8_17_hamming_reference_to_part_73():
    assert close(predicted("HAMMING_REFERENCE_TO_PART_73"), 6)


def test_8_18_hamming_flags():
    assert close(predicted("HAMMING_FLAGS"), 2)


def test_8_19_jaccard_recipe_winner():
    assert predicted("JACCARD_RECIPE_WINNER") == "Shortbread"


def test_8_20_cosine_recipe_winner():
    assert predicted("COSINE_RECIPE_WINNER") == "Sachertorte"


def test_8_21_covariance_of_readings():
    assert close(predicted("COVARIANCE_OF_READINGS"),
                 [[7.5, 7.0], [7.0, 7.5]])


def test_8_22_probes_are_equidistant_under_euclidean():
    assert predicted("PROBES_EQUIDISTANT_UNDER_EUCLIDEAN") is True


def test_8_23_mahalanobis_to_probe_across():
    assert close(predicted("MAHALANOBIS_TO_PROBE_ACROSS"), 6.0)


def test_8_24_raw_bearing_winner():
    assert predicted("RAW_BEARING_WINNER") == "R"


def test_8_25_standardised_bearing_winner():
    assert predicted("STANDARDISED_BEARING_WINNER") == "P"

#!/usr/bin/env bash
# Tests for the Day 107 lab. Run from the lab directory:
#   bash tests/run_tests.sh
#
# The harness proves the lesson's claims by running code and reading real
# values, never by reading source:
#
#   * one query, three candidates, and L1, L2 and cosine naming three
#     DIFFERENT winners -- the disagreement the whole day is built on;
#   * the p-norm of (3, 4) is 7, 5 and 4 at p = 1, 2 and infinity, falls
#     monotonically in between, and matches numpy.linalg.norm(v, ord=p) to
#     1e-12 for every p tested;
#   * the L1, L2 and L-infinity unit balls are strictly nested, and counting
#     grid cells inside the L2 one recovers pi;
#   * all four norm axioms hold for L1, L2 and L-infinity, and SQUARED
#     Euclidean distance fails absolute homogeneity by a factor of 2;
#   * the triangle inequality holds for L1, L2 and L-infinity in all six
#     orderings, and FAILS for cosine distance on 326 of 3375 binary triples,
#     with the two-dimensional counter-example named;
#   * Jaccard distance and Hamming distance satisfy it on all 4096 triples;
#   * Chebyshev accepts a part that L1 and L2 both rank as the worse one,
#     because the acceptance rule is an L-infinity ball;
#   * Jaccard and cosine rank the same two sets in OPPOSITE orders;
#   * two points the same Euclidean distance from the mean are 1.1142 and
#     6.0 apart under Mahalanobis, and the pure-Python Gauss-Jordan inverse
#     and numpy.linalg.inv disagree in the last bit while both round to 6;
#   * standardising two columns in mismatched units changes the winner, and a
#     UNIT change alone changes it too;
#   * nothing is downloaded, nothing is written outside the lab, and nothing
#     is left behind on disk.
#
# Everything runs offline. Nothing binds a port, nothing writes outside the
# lab or a temporary directory, nothing needs a key. Deterministic,
# non-interactive, exits 0 only if every check passes.
set -u

export PYTHONDONTWRITEBYTECODE=1

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

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

failures=0
checks=0

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

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

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

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

# The Python that owns that pytest is the one with numpy 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 107 — Choose Your Distance on Purpose"
echo

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

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

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

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

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

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

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

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

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

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

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

# The import guard. Both directories contain modules called `measures` and
# `catalogue`, 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 itertools
import math

import numpy as np

import catalogue
import measures

TOL = measures.TOL
Q = catalogue.QUERY
A = catalogue.ARTICLES

# -- the opening disagreement
print("l1_winner", measures.winner(Q, A, measures.l1_distance))
print("l2_winner", measures.winner(Q, A, measures.l2_distance))
print("linf_winner", measures.winner(Q, A, measures.linf_distance))
print("cosine_winner",
      measures.winner(Q, A, measures.cosine_similarity, higher_is_better=True))
print("distinct_winners", len({
    measures.winner(Q, A, measures.l1_distance),
    measures.winner(Q, A, measures.l2_distance),
    measures.winner(Q, A, measures.cosine_similarity, higher_is_better=True),
}))
print("l1_values", [measures.l1_distance(Q, v) for v in A.values()])
print("l2_aisle_exact", measures.l2_distance(Q, A["Aisle"]) == 5.0)
print("cosine_cartogram_is_one",
      abs(measures.cosine_similarity(Q, A["Cartogram"]) - 1.0) <= TOL)
print("cartogram_worst_l2",
      measures.rank(Q, A, measures.l2_distance)[-1][0])

# -- the p-norm family
v = (3.0, 4.0)
print("p_norms", [measures.p_norm(v, p) for p in (1, 2, math.inf)])
sweep = [measures.p_norm(v, p) for p in (1, 1.5, 2, 3, 4, 8, 16, 64)]
print("p_sweep_monotone", all(b <= a + TOL for a, b in zip(sweep, sweep[1:])))
worst_norm = max(
    abs(measures.p_norm(v, math.inf if np.isinf(p) else p)
        - float(np.linalg.norm(np.asarray(v), ord=p)))
    for p in (1, 1.5, 2, 3, 8, np.inf))
print("p_norm_vs_numpy", worst_norm <= TOL)
try:
    measures.p_norm(v, 0.5)
except ValueError:
    print("p_norm_refuses_below_one", True)
else:
    print("p_norm_refuses_below_one", False)

# -- the unit balls, counted on the same grid the script draws
W, H = 61, 25
counts = {1: 0, 2: 0, 3: 0}
for row in range(H):
    y = 1.25 - 2.5 * row / (H - 1)
    for col in range(W):
        x = -1.25 + 2.5 * col / (W - 1)
        if measures.p_norm((x, y), 1) <= 1.0:
            counts[1] += 1
        if measures.p_norm((x, y), 2) <= 1.0:
            counts[2] += 1
        if measures.p_norm((x, y), math.inf) <= 1.0:
            counts[3] += 1
print("ball_counts", counts[1], counts[2], counts[3])
print("balls_nested", counts[1] < counts[2] < counts[3])
cell = (2.5 / (W - 1)) * (2.5 / (H - 1))
print("l2_ball_area_recovers_pi", abs(counts[2] * cell - math.pi) < 0.25)

# -- the norm axioms
av, k = catalogue.AXIOM_VECTOR, catalogue.AXIOM_SCALAR
w = catalogue.TRIANGLE_TRIPLE[1]
homog, tri, zero_ok = True, True, True
for fn in (measures.l1_norm, measures.l2_norm, measures.linf_norm):
    homog &= abs(fn([k * x for x in av]) - abs(k) * fn(av)) <= TOL
    tri &= fn([a + b for a, b in zip(av, w)]) <= fn(av) + fn(w) + TOL
    zero_ok &= fn((0.0, 0.0, 0.0)) == 0.0 and fn(av) > 0.0
print("norm_homogeneity", homog)
print("norm_triangle", tri)
print("norm_zero_only_at_zero", zero_ok)
sq = sum(x * x for x in av)
print("squared_l2_doubles_to", sum((2 * x) ** 2 for x in av) / sq)

# -- metrics and non-metrics
d = measures.cosine_distance
east, diag, north = catalogue.EAST, catalogue.DIAGONAL, catalogue.NORTH
print("cosine_detour", round(d(east, diag) + d(diag, north), 6))
print("cosine_direct", round(d(east, north), 6))
print("cosine_violates", d(east, diag) + d(diag, north) < d(east, north) - TOL)
print("cosine_zero_between_different_vectors",
      abs(d((1.0, 0.0), (2.0, 0.0))) <= TOL)
ok = True
for fn in (measures.l1_distance, measures.l2_distance, measures.linf_distance):
    for a, b, c in itertools.permutations(catalogue.TRIANGLE_TRIPLE):
        ok &= fn(a, b) + fn(b, c) >= fn(a, c) - TOL
print("lp_triangle_all_orderings", ok)

universe = ("a", "b", "c", "d")
subsets = [frozenset(c) for r in range(len(universe) + 1)
           for c in itertools.combinations(universe, r)]
jac_ok = all(
    measures.jaccard_distance(a, b) + measures.jaccard_distance(b, c)
    >= measures.jaccard_distance(a, c) - TOL
    for a, b, c in itertools.product(subsets, repeat=3))
print("jaccard_triples", len(subsets) ** 3)
print("jaccard_is_metric", jac_ok)

bits = list(itertools.product((0, 1), repeat=4))
ham_ok = all(
    measures.hamming_distance(a, b) + measures.hamming_distance(b, c)
    >= measures.hamming_distance(a, c)
    for a, b, c in itertools.product(bits, repeat=3))
print("hamming_is_metric", ham_ok)

vectors = [x for x in itertools.product((0, 1), repeat=4) if any(x)]
violations = sum(
    1 for a, b, c in itertools.product(vectors, repeat=3)
    if d(a, b) + d(b, c) < d(a, c) - TOL)
print("cosine_violations", violations, len(vectors) ** 3)

# -- the four data shapes
fa, fb = catalogue.FLOOR_FROM, catalogue.FLOOR_TO
print("warehouse", measures.l1_distance(fa, fb), measures.l2_distance(fa, fb),
      measures.linf_distance(fa, fb))
nominal, limit = catalogue.NOMINAL_PART, catalogue.PART_TOLERANCE_MM
pa = catalogue.MEASURED_PARTS["batch-A"]
pb = catalogue.MEASURED_PARTS["batch-B"]
print("batch_a_passes", measures.linf_distance(pa, nominal) <= limit + TOL)
print("batch_b_fails", measures.linf_distance(pb, nominal) > limit)
print("l1_prefers_batch_b",
      measures.l1_distance(pb, nominal) < measures.l1_distance(pa, nominal))
print("l2_prefers_batch_b",
      measures.l2_distance(pb, nominal) > measures.l2_distance(pa, nominal))

ref = catalogue.REFERENCE_RECORD
print("hamming_records", [measures.hamming_distance(ref, r)
                          for r in catalogue.CANDIDATE_RECORDS.values()])
print("hamming_flags",
      measures.hamming_distance(catalogue.FLAGS_A, catalogue.FLAGS_B))
print("hamming_equals_l1_on_bits",
      abs(measures.l1_distance(catalogue.FLAGS_A, catalogue.FLAGS_B)
          - measures.hamming_distance(catalogue.FLAGS_A,
                                      catalogue.FLAGS_B)) <= TOL)

rq = catalogue.RECIPE_QUERY
axes = measures.vocabulary(rq, *catalogue.RECIPES.values())
qv = measures.to_binary_vector(rq, axes)
jac = {n: measures.jaccard_similarity(rq, s)
       for n, s in catalogue.RECIPES.items()}
cos = {n: measures.cosine_similarity(qv, measures.to_binary_vector(s, axes))
       for n, s in catalogue.RECIPES.items()}
print("jaccard_recipe_winner", max(jac, key=jac.get))
print("cosine_recipe_winner", max(cos, key=cos.get))
print("recipe_winners_differ", max(jac, key=jac.get) != max(cos, key=cos.get))
print("jaccard_values", round(jac["Sachertorte"], 6), round(jac["Shortbread"], 6))
print("cosine_values", round(cos["Sachertorte"], 6), round(cos["Shortbread"], 6))

# -- Mahalanobis
data = catalogue.SENSOR_READINGS
mean = measures.column_means(data)
cov = measures.covariance_matrix(data)
inv = measures.inverse(cov)
print("covariance", cov)
print("covariance_det", round(cov[0][0] * cov[1][1] - cov[0][1] * cov[1][0], 10))
print("inverse_matches_numpy",
      float(np.max(np.abs(np.asarray(inv) - np.linalg.inv(np.asarray(cov)))))
      <= TOL)
eu_a = measures.l2_distance(catalogue.PROBE_ALONG, mean)
eu_x = measures.l2_distance(catalogue.PROBE_ACROSS, mean)
print("probes_equidistant_euclidean", abs(eu_a - eu_x) <= TOL)
print("euclidean_probe_distance", round(eu_a, 6))
ma_a = measures.mahalanobis_distance(catalogue.PROBE_ALONG, mean, inv)
ma_x = measures.mahalanobis_distance(catalogue.PROBE_ACROSS, mean, inv)
print("mahalanobis_along", round(ma_a, 6))
print("mahalanobis_across", round(ma_x, 6))
print("mahalanobis_across_repr", repr(ma_x))
inv_np = np.linalg.inv(np.asarray(cov))
z = np.asarray(catalogue.PROBE_ACROSS)
print("mahalanobis_across_numpy_repr", repr(float(math.sqrt(z @ inv_np @ z))))
print("mahalanobis_routes_agree_within_tol",
      abs(ma_x - float(math.sqrt(z @ inv_np @ z))) <= TOL)
identity = [[1.0, 0.0], [0.0, 1.0]]
print("identity_gives_euclidean", max(
    abs(measures.mahalanobis_distance(p, mean, identity)
        - measures.l2_distance(p, mean))
    for p in (catalogue.PROBE_ALONG, catalogue.PROBE_ACROSS,
              (-2.5, 4.75))) <= TOL)
print("eigenvalues", sorted(round(float(x), 10)
                            for x in np.linalg.eigvalsh(np.asarray(cov))))

# -- scaling
B, BQ = catalogue.BEARINGS, catalogue.BEARING_QUERY
rows = list(B.values())
raw_order = [n for n, _ in measures.rank(BQ, B, measures.l2_distance)]
means, stds = measures.column_means(rows), measures.column_stds(rows)
qs = measures.standardise([BQ], means, stds)[0]
bs = {n: measures.standardise([v], means, stds)[0] for n, v in B.items()}
std_order = [n for n, _ in measures.rank(qs, bs, measures.l2_distance)]
print("raw_bearing_order", raw_order)
print("standardised_bearing_order", std_order)
print("bearing_winner_changed", raw_order[0] != std_order[0])
print("bore_share_max", "%.3e" % max(
    (BQ[0] - r[0]) ** 2 / ((BQ[0] - r[0]) ** 2 + (BQ[1] - r[1]) ** 2)
    for r in rows if r != tuple(BQ)))
micro_q = (BQ[0] * 1e6, BQ[1])
micro = {n: (v[0] * 1e6, v[1]) for n, v in B.items()}
print("unit_change_winner",
      measures.winner(micro_q, micro, measures.l2_distance))
maha_order = [n for n, _ in measures.rank(
    BQ, B, lambda a, b: measures.mahalanobis_distance(
        a, b, measures.inverse(measures.covariance_matrix(rows))))]
print("mahalanobis_bearing_order", maha_order)
print("cosine_not_unit_invariant",
      [n for n, _ in measures.rank(BQ, B, measures.cosine_similarity,
                                   higher_is_better=True)]
      != [n for n, _ in measures.rank(micro_q, micro,
                                      measures.cosine_similarity,
                                      higher_is_better=True)])
doubled = {n: tuple(2 * c for c in v) for n, v in B.items()}
print("cosine_vector_scale_invariant",
      [n for n, _ in measures.rank(BQ, B, measures.cosine_similarity,
                                   higher_is_better=True)]
      == [n for n, _ in measures.rank(BQ, doubled, measures.cosine_similarity,
                                      higher_is_better=True)])

# -- the seeded sweep
rng = np.random.default_rng(107)
spread = np.array([0.04, 500.0])
flips = 0
for _ in range(2000):
    cat = rng.random((6, 2)) * spread
    query = rng.random(2) * spread
    raw_best = int(np.argmin(np.linalg.norm(cat - query, axis=1)))
    mu, sd = cat.mean(axis=0), cat.std(axis=0)
    std_best = int(np.argmin(
        np.linalg.norm((cat - mu) / sd - (query - mu) / sd, axis=1)))
    flips += raw_best != std_best
print("seeded_sweep_flips", flips)
print("seeded_sweep_in_range", 0.35 <= flips / 2000 <= 0.75)

# -- guard rails
try:
    measures.l2_distance((1.0, 2.0, 3.0), (1.0, 2.0))
except measures.DimensionMismatch:
    print("length_mismatch_raises", True)
else:
    print("length_mismatch_raises", False)
try:
    measures.cosine_similarity((0.0, 0.0), (1.0, 1.0))
except ValueError:
    print("zero_vector_raises", True)
else:
    print("zero_vector_raises", False)
try:
    measures.inverse([[1.0, 2.0], [2.0, 4.0]])
except ValueError:
    print("singular_matrix_raises", True)
else:
    print("singular_matrix_raises", False)
print("ties_break_by_name", [n for n, _ in measures.rank(
    (0.0, 0.0), {"zulu": (1.0, 1.0), "alpha": (1.0, 1.0), "mike": (1.0, 1.0)},
    measures.l2_distance)])
PY
)"

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

check_eq "L1 picks Aisle" "Aisle" "$(get l1_winner)"
check_eq "L2 picks Beacon" "Beacon" "$(get l2_winner)"
check_eq "L-infinity also picks Beacon" "Beacon" "$(get linf_winner)"
check_eq "cosine picks Cartogram" "Cartogram" "$(get cosine_winner)"
check_eq "three measures name three DIFFERENT winners" "3" \
  "$(get distinct_winners)"
check_eq "the L1 distances are 5, 6 and 20" "[5, 6, 20]" "$(get l1_values)"
check_eq "the L2 distance to Aisle is exactly 5" "True" "$(get l2_aisle_exact)"
check_eq "the cosine to Cartogram is 1 within tolerance" "True" \
  "$(get cosine_cartogram_is_one)"
check_eq "and Cartogram is the WORST answer under L2" "Cartogram" \
  "$(get cartogram_worst_l2)"

check_eq "the p-norm of (3, 4) is 7, 5 and 4 at p = 1, 2 and infinity" \
  "[7.0, 5.0, 4.0]" "$(get p_norms)"
check_eq "the p-norm falls as p rises" "True" "$(get p_sweep_monotone)"
check_eq "and matches numpy.linalg.norm(v, ord=p) within 1e-12" "True" \
  "$(get p_norm_vs_numpy)"
check_eq "p below 1 is refused rather than answered" "True" \
  "$(get p_norm_refuses_below_one)"
check_eq "the three unit balls contain 469, 723 and 931 grid cells" \
  "469 723 931" "$(get ball_counts)"
check_eq "so the L1 ball sits inside L2 sits inside L-infinity" "True" \
  "$(get balls_nested)"
check_eq "and counting cells inside the L2 ball recovers pi" "True" \
  "$(get l2_ball_area_recovers_pi)"

check_eq "absolute homogeneity holds for L1, L2 and L-infinity" "True" \
  "$(get norm_homogeneity)"
check_eq "the triangle inequality holds for all three" "True" \
  "$(get norm_triangle)"
check_eq "and each is zero only at the zero vector" "True" \
  "$(get norm_zero_only_at_zero)"
check_eq "doubling a vector QUADRUPLES squared Euclidean distance" "4.0" \
  "$(get squared_l2_doubles_to)"

check_eq "cosine distance via the diagonal costs 0.585786" "0.585786" \
  "$(get cosine_detour)"
check_eq "and going direct costs 1.0, which is more" "1.0" \
  "$(get cosine_direct)"
# Section 6 re-runs this script with D107_SELF_TEST=1, which swaps ONE
# expectation below for a deliberately wrong one -- the naive belief that
# cosine distance behaves like a distance. That is how the harness proves it
# can fail rather than merely asserting that it could.
expected_cosine_violates="True"
if [ -n "${D107_SELF_TEST:-}" ]; then
  expected_cosine_violates="False"   # the naive belief, deliberately wrong here
fi
check_eq "so cosine distance VIOLATES the triangle inequality" \
  "${expected_cosine_violates}" "$(get cosine_violates)"
check_eq "cosine distance is also 0 between two different vectors" "True" \
  "$(get cosine_zero_between_different_vectors)"
check_eq "L1, L2 and L-infinity hold in all six orderings" "True" \
  "$(get lp_triangle_all_orderings)"
check_eq "Jaccard distance was checked on 4096 triples" "4096" \
  "$(get jaccard_triples)"
check_eq "and is a metric on every one of them" "True" \
  "$(get jaccard_is_metric)"
check_eq "Hamming distance is a metric on all 4096 bit triples" "True" \
  "$(get hamming_is_metric)"
check_eq "cosine distance fails on 326 of 3375 binary triples" "326 3375" \
  "$(get cosine_violations)"

check_eq "one displacement gives 14, 10 and 8" "14.0 10.0 8.0" \
  "$(get warehouse)"
check_eq "Chebyshev ACCEPTS batch-A" "True" "$(get batch_a_passes)"
check_eq "and REJECTS batch-B" "True" "$(get batch_b_fails)"
check_eq "even though L1 ranks batch-B as the better part" "True" \
  "$(get l1_prefers_batch_b)"
check_eq "and L2 ranks it the other way, so L1 and L2 disagree too" "True" \
  "$(get l2_prefers_batch_b)"

check_eq "Hamming on the parts register gives 1, 3 and 6" "[1, 3, 6]" \
  "$(get hamming_records)"
check_eq "and 2 on the bit flags" "2" "$(get hamming_flags)"
check_eq "on bits, Hamming equals L1 exactly" "True" \
  "$(get hamming_equals_l1_on_bits)"

check_eq "Jaccard prefers Shortbread" "Shortbread" \
  "$(get jaccard_recipe_winner)"
check_eq "cosine prefers Sachertorte on the SAME two sets" "Sachertorte" \
  "$(get cosine_recipe_winner)"
check_eq "so the two rank set data in opposite orders" "True" \
  "$(get recipe_winners_differ)"
check_eq "Jaccard: 4/11 against 2/5" "0.363636 0.4" "$(get jaccard_values)"
check_eq "cosine: 4/sqrt(44) against 2/sqrt(12)" "0.603023 0.57735" \
  "$(get cosine_values)"

check_eq "the covariance of the readings is exactly [[7.5, 7], [7, 7.5]]" \
  "[[7.5, 7.0], [7.0, 7.5]]" "$(get covariance)"
check_eq "its determinant is exactly 7.25" "7.25" "$(get covariance_det)"
check_eq "the pure-Python inverse matches numpy.linalg.inv" "True" \
  "$(get inverse_matches_numpy)"
check_eq "both probes are the same Euclidean distance from the mean" "True" \
  "$(get probes_equidistant_euclidean)"
check_eq "that distance is sqrt(18) = 4.242641" "4.242641" \
  "$(get euclidean_probe_distance)"
check_eq "Mahalanobis says 1.114172 ALONG the grain" "1.114172" \
  "$(get mahalanobis_along)"
check_eq "and 6.0 ACROSS it" "6.0" "$(get mahalanobis_across)"
check_eq "Gauss-Jordan gives exactly 6.0" "6.0" \
  "$(get mahalanobis_across_repr)"
check_eq "numpy.linalg.inv gives 5.999999999999999 for the same quantity" \
  "5.999999999999999" "$(get mahalanobis_across_numpy_repr)"
check_eq "the two routes agree within the stated tolerance" "True" \
  "$(get mahalanobis_routes_agree_within_tol)"
check_eq "substituting the identity gives back Euclidean exactly" "True" \
  "$(get identity_gives_euclidean)"
check_eq "the covariance eigenvalues are 0.5 and 14.5" "[0.5, 14.5]" \
  "$(get eigenvalues)"

check_eq "raw Euclidean ranks the bearings R, U, P, S, T, V" \
  "['R', 'U', 'P', 'S', 'T', 'V']" "$(get raw_bearing_order)"
check_eq "standardised, it ranks them P, U, R, S, T, V" \
  "['P', 'U', 'R', 'S', 'T', 'V']" "$(get standardised_bearing_order)"
check_eq "so standardising CHANGES the winner" "True" \
  "$(get bearing_winner_changed)"
check_eq "and before scaling the bore column contributes at most 0.0036%" \
  "3.600e-05" "$(get bore_share_max)"
check_eq "changing the bore unit ALONE also changes the winner" "P" \
  "$(get unit_change_winner)"
check_eq "Mahalanobis on the raw numbers gives a third answer" \
  "['U', 'P', 'S', 'T', 'R', 'V']" "$(get mahalanobis_bearing_order)"
check_eq "cosine is NOT invariant to a change of column units" "True" \
  "$(get cosine_not_unit_invariant)"
check_eq "cosine IS invariant to scaling a whole vector" "True" \
  "$(get cosine_vector_scale_invariant)"
check_eq "in 2000 seeded random catalogues the winner changed 1090 times" \
  "1090" "$(get seeded_sweep_flips)"
check_eq "which is between a third and three quarters" "True" \
  "$(get seeded_sweep_in_range)"

check_eq "comparing different lengths raises rather than truncating" "True" \
  "$(get length_mismatch_raises)"
check_eq "cosine of the zero vector raises rather than returning 0" "True" \
  "$(get zero_vector_raises)"
check_eq "a singular covariance refuses to invert" "True" \
  "$(get singular_matrix_raises)"
check_eq "ranking ties break by name, so runs are deterministic" \
  "['alpha', 'mike', 'zulu']" "$(get ties_break_by_name)"

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

# A green test suite proves nothing until you have watched it go red. This
# section re-runs the whole script with one expectation deliberately swapped
# for the naive belief that cosine distance satisfies the triangle inequality,
# and asserts that the re-run reports the failure and exits non-zero. If this
# section passes, section 5 is not decorative.
if [ -z "${D107_SELF_TEST:-}" ]; then
  self_out="$(D107_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: so cosine distance VIOLATES the triangle inequality"*)
      check "the failing check is named in the output with both values" "yes" ;;
    *) check "the failing check is named in the output with both values" "no" ;;
  esac
  case "${self_out}" in
    *", 1 failure(s)."*)
      check "the summary line counts exactly one failure" "yes" ;;
    *) check "the summary line counts exactly one failure" "no" ;;
  esac
else
  echo "  (self-test run: section 6 does not recurse)"
fi

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

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

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

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

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

# measures.py must not import NumPy. The whole evidential value of agreeing
# with numpy.linalg.norm depends on it, so it is checked rather than trusted.
if grep -qE '^\s*(import|from)\s+numpy' "${lab_dir}/examples/measures.py" \
     "${lab_dir}/starter/measures.py" 2>/dev/null; then
  check "measures.py computes without NumPy, so agreeing with it means something" "no"
else
  check "measures.py computes without NumPy, so agreeing with it means something" "yes"
fi

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