Hands-on lab — Day 109: Partial Derivatives and Gradients

"""What a partial derivative is: one input moves, the rest are held still.

Run from inside `examples/`:

    ../.venv/bin/python3 01_hold_everything_else_still.py

Every claim printed here is asserted. If the script exits 0 and ends with
"every assertion held.", the numbers above it were computed, not typed.
"""

from __future__ import annotations

import numpy as np

import surfaces as S
from gradients import partial

print(__doc__.splitlines()[0])
print()

# --------------------------------------------------------------------------
print("1. A function of two inputs, and one slice through it")
# --------------------------------------------------------------------------
#
# f(x, y) = x^2 + 3y^2. Stand at (2, 1). There is no single "the slope" here,
# because there is no single direction to walk in. There is a slope if you
# walk east, a different slope if you walk north, and a different one again
# for every bearing in between.
#
# A partial derivative picks one of those and only one: freeze every input but
# the chosen one, which turns the function of two variables into a function of
# ONE variable, and then take Day 108's ordinary derivative of that.

point = (2.0, 1.0)
print(f"  f(x, y) = x^2 + 3y^2, standing at (x, y) = {point}")
print(f"  f at that point                    {S.bowl(point):.6f}")
print()
print("  Freeze y at 1. What is left is a function of x alone:")
print("      g(x) = f(x, 1) = x^2 + 3")
for x in (1.0, 1.5, 2.0, 2.5, 3.0):
    print(f"        g({x:>3}) = {S.bowl((x, 1.0)):8.4f}")
print("  and dg/dx = 2x, which at x = 2 is 4.")
print()
print("  Now freeze x at 2 instead. What is left is a function of y alone:")
print("      k(y) = f(2, y) = 4 + 3y^2")
for y in (0.0, 0.5, 1.0, 1.5, 2.0):
    print(f"        k({y:>3}) = {S.bowl((2.0, y)):8.4f}")
print("  and dk/dy = 6y, which at y = 1 is 6.")
print()
print("  Those two numbers, 4 and 6, are the two partial derivatives at (2, 1).")
print("  They are written  df/dx = 4  and  df/dy = 6, with the rounded d")
print("  rather than the straight one, and the rounded d is the entire notice")
print("  that other inputs exist and are being held still.")

assert S.bowl(point) == 7.0
assert S.bowl_gradient(point)[0] == 4.0
assert S.bowl_gradient(point)[1] == 6.0

# --------------------------------------------------------------------------
print()
print("2. The same two numbers, measured instead of derived")
# --------------------------------------------------------------------------
#
# Nothing above needed a computer. But a computer cannot read x^2 + 3y^2 and
# differentiate it; it can only evaluate it. So do what Day 108 did: nudge the
# input a little each way and divide the change by the distance moved. The
# only new instruction is "and change nothing else".

h = S.H_DEFAULT
print(f"  step size h = {h:g}")
print()
for index, name in ((0, "x"), (1, "y")):
    base = np.asarray(point, dtype=float)
    up = base.copy()
    down = base.copy()
    up[index] += h
    down[index] -= h
    measured = partial(S.bowl, point, index, h)
    exact = S.bowl_gradient(point)[index]
    print(f"  df/d{name}:")
    shown_up = "(" + ", ".join(f"{v:.5f}" for v in up) + ")"
    shown_down = "(" + ", ".join(f"{v:.5f}" for v in down) + ")"
    print(f"    point nudged up      {shown_up}   f = {S.bowl(up):.12f}")
    print(f"    point nudged down    {shown_down}   f = {S.bowl(down):.12f}")
    print(f"    difference / (2h)    {measured:.12f}")
    print(f"    exact, by hand       {exact:.12f}")
    print(f"    error                {abs(measured - exact):.3e}")
    assert abs(measured - exact) < S.GRADIENT_TOL

print()
print("  Both errors are far below the tolerance this lab asserts,")
print(f"  which is {S.GRADIENT_TOL:g}. Section 3 explains why they are THIS small.")

# --------------------------------------------------------------------------
print()
print("3. Why the error is roundoff rather than method error, here")
# --------------------------------------------------------------------------
#
# A central difference has a truncation error proportional to h^2 times the
# third derivative. x^2 + 3y^2 has no third derivative worth the name -- it is
# identically zero -- so the h^2 term vanishes and what is left is only the
# floating-point noise from subtracting two nearly equal numbers.
#
# The algebra is short enough to show. For g(x) = x^2:
#
#     ((x+h)^2 - (x-h)^2) / (2h)
#   = (x^2 + 2xh + h^2 - x^2 + 2xh - h^2) / (2h)
#   = 4xh / (2h)
#   = 2x,    for ANY h, exactly.
#
# So for a quadratic the central difference is not an approximation at all. It
# is the answer, and the only thing between you and it is float64.

print("  Central difference on a quadratic is algebraically EXACT:")
print("    ((x+h)^2 - (x-h)^2) / (2h) = 4xh / (2h) = 2x, for any h at all.")
print()
print("  So changing h should barely move the answer. It does not:")
print(f"    {'h':>10}  {'df/dx at (2, 1)':>20}  {'error':>12}")
for k in (1, 2, 3, 4, 5, 6):
    hh = 10.0 ** (-k)
    value = partial(S.bowl, point, 0, hh)
    print(f"    {hh:10.0e}  {value:20.14f}  {abs(value - 4.0):12.3e}")
    assert abs(value - 4.0) < 1e-7

print()
print("  Script 06 does this on a genuine cubic, where the h^2 term is real,")
print("  and the same table becomes the U-shaped curve from Day 108.")

# --------------------------------------------------------------------------
print()
print("4. The function whose partials need the other variable")
# --------------------------------------------------------------------------
#
# x^2 + 3y^2 is a soft case: freezing y leaves a function of x with no y in it
# at all, so it is easy to believe the two variables were never really
# interacting. f(x, y) = xy destroys that comfort.

print("  f(x, y) = xy.  df/dx = y  and  df/dy = x.")
print("  The slope in x depends on where you are in Y. Walk along the x-axis,")
print("  where y = 0, and f is identically zero, so the slope in x is zero.")
print("  Step off that line and it stops being zero.")
print()
print(f"    {'point':>14}  {'df/dx exact':>12}  {'measured':>16}  {'df/dy exact':>12}  {'measured':>16}")
for p in ((1.0, 0.0), (1.0, 1.0), (1.0, 5.0), (3.0, -2.0)):
    ex = S.product_gradient(p)
    mx = partial(S.product, p, 0)
    my = partial(S.product, p, 1)
    print(f"    {str(p):>14}  {ex[0]:12.4f}  {mx:16.12f}  {ex[1]:12.4f}  {my:16.12f}")
    assert abs(mx - ex[0]) < S.GRADIENT_TOL
    assert abs(my - ex[1]) < S.GRADIENT_TOL

print()
print("  Read the first row again: at (1, 0) the slope in x is exactly zero,")
print("  and the surface is emphatically not flat there -- the slope in y is 1.")
print("  A single partial derivative being zero says nothing about the point.")

# --------------------------------------------------------------------------
print()
print("5. Every input gets one, however many there are")
# --------------------------------------------------------------------------
#
# Nothing above used the fact that there were two inputs. `partial` takes an
# index, so it works on a point of any length. Here is a function of three.

three = S.START_PARAMS
print(f"  A three-parameter loss at {three}: L = {S.model_loss(three)}")
exact3 = S.model_loss_gradient(three)
for i, label in enumerate(("w1", "w2", "c")):
    measured = partial(S.model_loss, three, i)
    print(f"    dL/d{label:<2}  exact {exact3[i]:8.4f}   measured {measured:18.12f}"
          f"   error {abs(measured - exact3[i]):.3e}")
    assert abs(measured - exact3[i]) < S.GRADIENT_TOL

print()
print("  Three inputs, three partial derivatives, and six evaluations of L to")
print("  get them -- two per input. Script 07 follows that cost to its")
print("  conclusion, which is the reason autodiff exists.")

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

"""Collect every partial derivative into one vector and you have the gradient.

Run from inside `examples/`:

    ../.venv/bin/python3 02_the_gradient_vector.py
"""

from __future__ import annotations

import numpy as np

import surfaces as S
from gradients import angle_degrees, gradient, magnitude, unit

print(__doc__.splitlines()[0])
print()

# --------------------------------------------------------------------------
print("1. Six surfaces, five points each, numerical against exact")
# --------------------------------------------------------------------------
#
# The gradient is not a new idea. It is the partial derivatives of script 01
# written side by side in square brackets instead of on separate lines. What
# makes it worth a name is that the result is a VECTOR, and Days 99 to 103
# already taught what to do with one: it has a length, it has a direction, and
# it can be dotted with another vector.
#
# Everything below is checked against a gradient worked out with a pencil.

worst = 0.0
count = 0
for name, (f, exact_gradient, expression, gradient_expression) in S.SURFACES.items():
    print(f"  f(x, y) = {expression:<14}   {gradient_expression}")
    print(f"    {'point':>14}  {'numerical gradient':>34}  {'exact':>18}  {'max error':>11}")
    for p in S.PROBE_POINTS:
        numeric = gradient(f, p)
        exact = exact_gradient(p)
        error = float(np.max(np.abs(numeric - exact)))
        worst = max(worst, error)
        count += 1
        shown_numeric = "[" + ", ".join(f"{v:14.10f}" for v in numeric) + "]"
        shown_exact = "[" + ", ".join(f"{v:7.3f}" for v in exact) + "]"
        print(f"    {str(p):>14}  {shown_numeric:>34}  {shown_exact:>18}  {error:11.3e}")
        assert error < S.GRADIENT_TOL, (name, p, error)
    print()

print(f"  {count} gradients checked. Worst single error {worst:.3e},")
print(f"  against an asserted tolerance of {S.GRADIENT_TOL:g}.")

# --------------------------------------------------------------------------
print()
print("2. The gradient is a vector, so it has a length and a bearing")
# --------------------------------------------------------------------------
#
# Two numbers come out of the same object and they answer different questions.
# The DIRECTION answers "which way is uphill". The LENGTH answers "how steep is
# it that way", in units of f gained per unit of distance walked.

print(f"  {'surface':>9}  {'point':>13}  {'gradient':>17}  {'length':>10}  {'bearing':>9}")
for name in ("bowl", "plane", "saddle", "cubic"):
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    for p in ((1.0, 1.0), (2.0, -1.0)):
        g = gradient(f, p)
        shown = "[" + ", ".join(f"{v:7.4f}" for v in g) + "]"
        print(f"  {name:>9}  {str(p):>13}  {shown:>17}  {magnitude(g):10.6f}"
              f"  {angle_degrees(g):8.3f}d")
        assert abs(magnitude(g) - magnitude(exact_gradient(p))) < S.GRADIENT_TOL

print()
print("  Read the bowl's two rows. At (1, 1) the gradient is about [2, 6]:")
print("  three times as much climb per step north as per step east, because")
print("  the 3 in front of y^2 makes the bowl three times steeper that way.")

# --------------------------------------------------------------------------
print()
print("3. A gradient is a direction in the INPUT space, not a point on the surface")
# --------------------------------------------------------------------------
#
# This is the most common way to misread the object. f(x, y) = x^2 + 3y^2 has
# a two-dimensional input and a one-dimensional output, and lives naturally as
# a surface in three dimensions. The gradient has TWO components, not three.
# It is an arrow drawn on the flat map you are standing on, not an arrow
# pointing up out of the hillside.

p = (1.0, 1.0)
g = gradient(S.bowl, p)
print(f"  f takes a point with {len(p)} coordinates and returns 1 number:")
print(f"    f{p} = {S.bowl(p)}")
print(f"  and its gradient has {g.size} components, matching the INPUT:")
print(f"    grad f{p} = [{g[0]:.6f}, {g[1]:.6f}]")
stepped = np.array(p) + g
print("  Add the gradient to the point and you get another point, "
      f"({stepped[0]:.4f}, {stepped[1]:.4f}),")
print("  which is a legal thing to do and is exactly what Day 111 will do,")
print("  with a minus sign in front.")
assert g.size == len(p)

# --------------------------------------------------------------------------
print()
print("4. The unit gradient: direction with the steepness divided out")
# --------------------------------------------------------------------------

for p in ((1.0, 1.0), (0.25, 0.75), (3.0, 0.5)):
    g = gradient(S.bowl, p)
    u = unit(g)
    print(f"  at {str(p):>13}  gradient [{g[0]:9.5f}, {g[1]:9.5f}]"
          f"   unit [{u[0]:8.5f}, {u[1]:8.5f}]   |unit| = {magnitude(u):.15f}")
    assert abs(magnitude(u) - 1.0) < 1e-12
    assert abs(angle_degrees(u) - angle_degrees(g)) < 1e-9

print()
print("  Same bearing, length exactly 1. Script 03 needs the unit version,")
print("  because a rate of change 'in a direction' is meaningless until the")
print("  direction has a fixed length -- otherwise drawing a longer arrow")
print("  would make the hill steeper.")

# --------------------------------------------------------------------------
print()
print("5. The zero gradient, and what it does not tell you")
# --------------------------------------------------------------------------

print("  Three different surfaces, all with gradient [0, 0] at the origin:")
print(f"    {'surface':>9}  {'gradient at origin':>32}  {'length':>12}  what the origin IS")
for name, kind, why in S.STATIONARY_AT_ORIGIN:
    f = S.SURFACES[name][0]
    g = gradient(f, (0.0, 0.0))
    shown = "[" + ", ".join(f"{v:14.11f}" for v in g) + "]"
    print(f"    {name:>9}  {shown:>32}  {magnitude(g):12.3e}  {kind} -- {why}")
    assert magnitude(g) < S.GRADIENT_TOL

print()
print("  The gradient is identical in all three cases and the points are not")
print("  remotely alike. A zero gradient says 'the ground is level here'. It")
print("  does not say whether you are at the bottom of a valley, on top of a")
print("  hill, or in a mountain pass. Script 05 shows the difference by")
print("  walking away from each one.")

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

"""The gradient really is the steepest way up. Measured, not asserted.

Run from inside `examples/`:

    ../.venv/bin/python3 03_steepest_ascent.py
"""

from __future__ import annotations

import numpy as np

import surfaces as S
from gradients import (
    angle_degrees,
    angular_gap_degrees,
    directional_derivative,
    directional_derivative_direct,
    gradient,
    magnitude,
    sweep_directions,
    unit,
)

print(__doc__.splitlines()[0])
print()

# --------------------------------------------------------------------------
print("1. A directional derivative, measured two ways that must agree")
# --------------------------------------------------------------------------
#
# "How fast does f change if I walk THIS way" is a question the partial
# derivatives do not directly answer, because they only know about the axes.
# There are two ways to answer it.
#
#   Directly: step h forward along the direction and h back along it, and
#   divide by 2h. That is Day 108's central difference with the step taken
#   along a diagonal rather than along an axis. It never forms a gradient.
#
#   Via the gradient: dot the gradient with the unit direction. That is Day
#   103's dot product, and it is not obvious that it should work.
#
# They agree. That agreement is the reason the gradient is worth assembling:
# two numbers, computed once, answer the question for every direction at once.

point = (1.0, 1.0)
g = gradient(S.bowl, point)
print(f"  f(x, y) = x^2 + 3y^2 at {point}, gradient about [{g[0]:.4f}, {g[1]:.4f}]")
print()
print(f"    {'direction':>16}  {'unit direction':>20}  {'via gradient':>15}  {'measured direct':>16}  {'gap':>10}")
directions = (
    (1.0, 0.0),
    (0.0, 1.0),
    (1.0, 1.0),
    (-1.0, 2.0),
    (3.0, -1.0),
    (-2.0, -5.0),
    (7.0, 0.5),
)
for d in directions:
    u = unit(np.array(d))
    via = directional_derivative(S.bowl, point, d)
    direct = directional_derivative_direct(S.bowl, point, d)
    shown_u = "[" + ", ".join(f"{v:8.5f}" for v in u) + "]"
    print(f"    {str(d):>16}  {shown_u:>20}  {via:15.9f}  {direct:16.9f}  {abs(via - direct):10.2e}")
    assert abs(via - direct) < S.GRADIENT_TOL

print()
print("  Look at the first two rows. Walking due east gives 2.0 and walking")
print("  due north gives 6.0 -- which are exactly the two partial derivatives.")
print("  A partial derivative is just the directional derivative along an axis.")
assert abs(directional_derivative(S.bowl, point, (1.0, 0.0)) - 2.0) < S.GRADIENT_TOL
assert abs(directional_derivative(S.bowl, point, (0.0, 1.0)) - 6.0) < S.GRADIENT_TOL

print()
print("  And note the last row: [7, 0.5] is a long arrow and [1, 0] is a short")
print("  one, but the answer depends only on the bearing, because both were")
print("  scaled to length 1 first.")
print()
print("  One row is worth stopping on. The direction (3, -1) gives a rate of")
print("  exactly zero, and it is not a coincidence: 3 times 2 plus -1 times 6")
print("  is 0, so that direction is perpendicular to the gradient. Walk that")
print("  way and, to first order, f does not change at all. Script 04 is about")
print("  what that means geometrically.")
assert abs(directional_derivative(S.bowl, point, (3.0, -1.0))) < S.GRADIENT_TOL

# --------------------------------------------------------------------------
print()
print("2. Try every direction and see which one wins")
# --------------------------------------------------------------------------
#
# The claim is that no direction climbs faster than the gradient's. So try
# them: 360 bearings evenly spaced around the circle, each one measured
# DIRECTLY with a central difference along that bearing, so the gradient is
# nowhere involved in producing the numbers being compared.

print(f"  {S.N_DIRECTIONS} directions, one per degree, measured directly.")
print()
print(f"  {'surface':>8}  {'point':>13}  {'best bearing':>13}  {'gradient bearing':>17}"
      f"  {'gap':>8}  {'best rate':>12}  {'|gradient|':>12}")

trials = (
    ("bowl", (1.0, 1.0)),
    ("bowl", (0.25, 0.75)),
    ("bowl", (3.0, -2.0)),
    ("product", (2.0, -1.0)),
    ("saddle", (1.5, 0.5)),
    ("cubic", (1.0, 1.0)),
    ("plane", (-2.0, 4.0)),
)
worst_gap = 0.0
for name, p in trials:
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    angles, rates = sweep_directions(f, p)
    best = int(np.argmax(rates))
    best_bearing = float(np.degrees(angles[best]))
    gradient_bearing = angle_degrees(exact_gradient(p))
    gap = angular_gap_degrees(best_bearing, gradient_bearing)
    worst_gap = max(worst_gap, gap)
    steepness = magnitude(exact_gradient(p))
    print(f"  {name:>8}  {str(p):>13}  {best_bearing:12.1f}d  {gradient_bearing:16.4f}d"
          f"  {gap:7.4f}d  {rates[best]:12.7f}  {steepness:12.7f}")
    assert gap <= S.ANGLE_TOL_DEGREES, (name, p, gap)
    # No direction may beat the gradient's own magnitude, ever.
    assert rates[best] <= steepness + S.GRADIENT_TOL

print()
print(f"  Worst gap across all {len(trials)} trials: {worst_gap:.4f} degrees,")
print(f"  against an asserted tolerance of {S.ANGLE_TOL_DEGREES} degree.")
print()
print("  The gap is not zero and cannot be. With one sample per degree, the")
print("  nearest sampled bearing to the true one is at most half a degree")
print("  away. The tolerance is that bound plus a little slack -- it is a")
print("  property of the sampling, not of the calculus.")

# --------------------------------------------------------------------------
print()
print("3. The winning rate is the gradient's LENGTH, and the shortfall is exactly cos(gap)")
# --------------------------------------------------------------------------
#
# Since the directional derivative is grad . u and both grad and u have fixed
# lengths, Day 103's geometric reading of the dot product applies unchanged:
#
#     grad . u = |grad| * |u| * cos(angle) = |grad| * cos(angle)
#
# So the best possible rate is |grad|, achieved when the angle is zero, and
# any other bearing gets |grad| times the cosine of how far off it is. That is
# not a rough statement; it is checkable to nine decimal places.

print(f"  {'surface':>8}  {'point':>13}  {'best rate / |grad|':>20}  {'cos(gap)':>18}  {'difference':>12}")
for name, p in trials:
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    angles, rates = sweep_directions(f, p)
    best = int(np.argmax(rates))
    gap = angular_gap_degrees(float(np.degrees(angles[best])),
                             angle_degrees(exact_gradient(p)))
    ratio = rates[best] / magnitude(exact_gradient(p))
    predicted = float(np.cos(np.radians(gap)))
    print(f"  {name:>8}  {str(p):>13}  {ratio:20.12f}  {predicted:18.12f}"
          f"  {abs(ratio - predicted):12.2e}")
    assert abs(ratio - predicted) < 1e-9

print()
print("  Five of the seven gaps above are the identical 0.4349 degrees, which")
print("  looks suspicious and is not. Those five gradients have bearings whose")
print("  fractional part is the same -- 26.5651, 71.5651, 116.5651, 296.5651,")
print("  341.5651 -- because they are all arctangents of ratios of the same")
print("  small whole numbers, separated by exact multiples of 45 degrees. A")
print("  grid sampled every whole degree therefore misses each of them by the")
print("  same amount. The two rows that break the pattern, at bearings 83.6598")
print("  and 326.3099, have different gaps. Pick a point with less tidy")
print("  coordinates and the gap changes again:")
for p in ((1.0, 0.4), (2.3, 1.7)):
    angles, rates = sweep_directions(S.bowl, p)
    best = int(np.argmax(rates))
    gap = angular_gap_degrees(float(np.degrees(angles[best])),
                             angle_degrees(S.bowl_gradient(p)))
    print(f"    bowl at {str(p):>12}: gradient bearing"
          f" {angle_degrees(S.bowl_gradient(p)):8.4f}d, gap {gap:.4f}d")
    assert gap <= S.ANGLE_TOL_DEGREES

# --------------------------------------------------------------------------
print()
print("4. The other end: the worst direction, and the two that go nowhere")
# --------------------------------------------------------------------------

p = (1.0, 1.0)
angles, rates = sweep_directions(S.bowl, p)
best = int(np.argmax(rates))
worst = int(np.argmin(rates))
flat = np.argsort(np.abs(rates))[:2]
steepness = magnitude(S.bowl_gradient(p))

print(f"  f(x, y) = x^2 + 3y^2 at {p}, |gradient| = {steepness:.7f}")
print(f"    steepest UP        bearing {np.degrees(angles[best]):6.1f}d"
      f"   rate {rates[best]:+12.7f}")
print(f"    steepest DOWN      bearing {np.degrees(angles[worst]):6.1f}d"
      f"   rate {rates[worst]:+12.7f}")
for i in sorted(flat):
    print(f"    no change at all   bearing {np.degrees(angles[i]):6.1f}d"
          f"   rate {rates[i]:+12.7f}")

print()
print("  The steepest descent is the exact opposite bearing, 180 degrees round,")
print("  and its rate is the negative of the steepest ascent. That symmetry is")
print("  the whole of Day 111 in one line: to go DOWN, step against the")
print("  gradient.")
opposite = angular_gap_degrees(float(np.degrees(angles[best])),
                               float(np.degrees(angles[worst])))
print(f"    angular separation of the two extremes: {opposite:.4f} degrees")
assert abs(opposite - 180.0) < 1e-9
assert abs(rates[best] + rates[worst]) < 1e-6

print()
print("  The two bearings where the rate is nearest zero are 90 degrees from")
print("  the gradient, one each way. Walking along either one keeps f the same")
print("  to first order -- which means they run along the contour. Script 04")
print("  makes that precise without going anywhere near this sweep.")
for i in flat:
    off = angular_gap_degrees(float(np.degrees(angles[i])),
                              angle_degrees(S.bowl_gradient(p)))
    assert abs(off - 90.0) <= S.ANGLE_TOL_DEGREES, off

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

"""The gradient is perpendicular to the contour through the point.

This is the geometric fact that makes every optimisation picture in the rest
of the course readable, so it is demonstrated rather than asserted -- and
demonstrated CAREFULLY, because the obvious way to demonstrate it is circular.

Run from inside `examples/`:

    ../.venv/bin/python3 04_perpendicular_to_the_contour.py
"""

from __future__ import annotations

import numpy as np

import surfaces as S
from gradients import contour_chord, gradient, magnitude, unit

print(__doc__.splitlines()[0])
print()

# --------------------------------------------------------------------------
print("1. The trap this script is built to avoid")
# --------------------------------------------------------------------------
#
# The lazy demonstration is: take the gradient, rotate it 90 degrees, call
# that "the contour direction", and observe that it is perpendicular to the
# gradient. That proves nothing whatever -- it is perpendicular because it was
# constructed to be.
#
# So every contour used below is an exact algebraic curve, derived on paper
# from the function alone, with the gradient nowhere in its derivation. And
# before any of it is used, the script CHECKS that f really is constant along
# each curve, so the reader is not asked to trust the algebra either.

print("  Three exactly parametrised contours, none derived from a gradient:")
print()
print("    bowl     x^2 + 3y^2 = L    x = sqrt(L) cos t,  y = sqrt(L/3) sin t")
print("    product  xy = L           x = t,  y = L/t")
print("    dome     -(x^2+y^2) = L   x = sqrt(-L) cos t,  y = sqrt(-L) sin t")
print()
print("  Substituting the first into x^2 + 3y^2 gives L cos^2 t + 3 (L/3) sin^2 t,")
print("  which is L (cos^2 t + sin^2 t) = L for every t. No gradient anywhere.")
print()
print("  Checked numerically at eight parameter values per curve:")
print(f"    {'surface':>9}  {'level L':>9}  {'max |f(point) - L| over the curve':>36}")
for name, (f, contour, level, _t0) in S.CONTOURS.items():
    ts = np.linspace(0.3, 2.6, 8)
    drift = max(abs(f(contour(level, t)) - level) for t in ts)
    print(f"    {name:>9}  {level:9.2f}  {drift:36.3e}")
    assert drift < 1e-12, (name, drift)

# --------------------------------------------------------------------------
print()
print("2. Step along the contour, and dot with the gradient")
# --------------------------------------------------------------------------
#
# Take a point p on the curve at parameter t, and a second point q at
# parameter t + delta. The unit vector from p to q is a CHORD of the contour.
# Dot it with the unit gradient at p.
#
# The result is not exactly zero and should not be expected to be, because a
# chord is not a tangent: it is tilted away from the tangent by an angle of
# roughly delta. So the honest evidence is not one small number -- it is the
# number shrinking in step with delta.

print(f"  {'surface':>9}  {'delta':>9}  {'unit gradient . unit chord':>28}  {'ratio to previous':>18}")
for name, (f, contour, level, t0) in S.CONTOURS.items():
    previous = None
    for k in (2, 3, 4, 5, 6):
        delta = 10.0 ** (-k)
        chord, p, _q, f_p, f_q = contour_chord(f, contour, level, t0, delta)
        g = unit(gradient(f, p))
        dot = float(np.dot(g, chord))
        ratio = "" if previous is None else f"{previous / abs(dot):18.4f}"
        print(f"  {name:>9}  {delta:9.0e}  {dot:+28.10e}  {ratio:>18}")
        assert abs(f_p - f_q) < 1e-12, "the two points are not on the same contour"
        if previous is not None:
            # Tenfold smaller step, tenfold smaller dot product: first order.
            assert 9.0 < previous / abs(dot) < 11.0, (name, delta, previous / abs(dot))
        previous = abs(dot)
    print()

print("  Every ratio is close to 10. Divide the step by ten and the dot")
print("  product divides by ten. That is what 'it goes to zero' looks like")
print("  when you can only ever take a finite step: not a small number, but a")
print("  number that shrinks at exactly the rate the geometry predicts.")

# --------------------------------------------------------------------------
print()
print("3. The same claim at the tolerance the lab actually asserts")
# --------------------------------------------------------------------------

print(f"  step along the contour: delta = {S.CONTOUR_DELTA:g}")
print(f"  asserted tolerance:            {S.CONTOUR_DOT_TOL:g}")
print()
print(f"  {'surface':>9}  {'t':>6}  {'point on the contour':>26}  {'unit gradient':>24}  {'dot':>13}")
worst = 0.0
for name, (f, contour, level, _t0) in S.CONTOURS.items():
    for t in (0.4, 0.9, 1.4, 1.9):
        chord, p, _q, _fp, _fq = contour_chord(f, contour, level, t, S.CONTOUR_DELTA)
        g = unit(gradient(f, p))
        dot = abs(float(np.dot(g, chord)))
        worst = max(worst, dot)
        shown_p = "(" + ", ".join(f"{v:10.6f}" for v in p) + ")"
        shown_g = "[" + ", ".join(f"{v:9.6f}" for v in g) + "]"
        print(f"  {name:>9}  {t:6.1f}  {shown_p:>26}  {shown_g:>24}  {dot:13.3e}")
        assert dot < S.CONTOUR_DOT_TOL, (name, t, dot)

print()
print(f"  Worst dot product across all {3 * 4} checks: {worst:.3e}, which is")
print(f"  {S.CONTOUR_DOT_TOL / worst:.0f} times inside the asserted tolerance. The tolerance was")
print("  chosen from section 2's measured rate before any of these ran, not")
print("  tightened afterwards until it looked impressive.")

# --------------------------------------------------------------------------
print()
print("4. Exactly zero, if you use the exact tangent instead of a chord")
# --------------------------------------------------------------------------
#
# The residual above is entirely the chord's fault. Differentiate the bowl's
# parametrisation with respect to t and you get the true tangent:
#
#     p(t)  = ( sqrt(L) cos t,       sqrt(L/3) sin t )
#     p'(t) = ( -sqrt(L) sin t,      sqrt(L/3) cos t )
#
# Dot that with the exact gradient (2x, 6y) = (2 sqrt(L) cos t, 6 sqrt(L/3) sin t):
#
#     -2L sin t cos t  +  6 (L/3) sin t cos t  =  (-2L + 2L) sin t cos t  =  0
#
# Identically zero, for every t and every level. No tolerance required.

level = 4.0
print(f"  bowl, contour level L = {level}")
print(f"    {'t':>6}  {'exact tangent':>26}  {'exact gradient':>26}  {'dot':>13}")
a = np.sqrt(level)
b = np.sqrt(level / 3.0)
for t in (0.0, 0.4, 0.9, 1.4, 1.9, 2.7):
    p = S.bowl_contour(level, t)
    tangent = np.array([-a * np.sin(t), b * np.cos(t)])
    g = S.bowl_gradient(p)
    dot = float(np.dot(tangent, g))
    shown_t = "[" + ", ".join(f"{v:10.6f}" for v in tangent) + "]"
    shown_g = "[" + ", ".join(f"{v:10.6f}" for v in g) + "]"
    print(f"    {t:6.1f}  {shown_t:>26}  {shown_g:>26}  {dot:13.3e}")
    assert abs(dot) < 1e-14

print()
print("  So the perpendicularity is exact and the small numbers in sections 2")
print("  and 3 are an artefact of measuring with finite steps, not a hedge.")

# --------------------------------------------------------------------------
print()
print("5. Why this is the fact that makes gradient descent make sense")
# --------------------------------------------------------------------------
#
# A contour is the set of points where f has one particular value -- a level
# set. On a real map it is the line joining points of equal height. Walking
# along it, you gain nothing and lose nothing.
#
# Perpendicular to it is therefore the only direction with anything to gain,
# and the gradient points along it. Everything else follows: the steepest
# ascent of script 03 is perpendicular to the contour, the steepest descent is
# the same line the other way, and the picture of an optimiser crossing
# contour lines at right angles -- which is what Day 112 will draw -- is not a
# stylisation. It is what the arrows do.

p = np.array([1.0, 1.0])
g = gradient(S.bowl, p)
along = np.array([-g[1], g[0]])   # a right angle to the gradient, by rotation
step = 0.001
print(f"  Stand at ({p[0]:.1f}, {p[1]:.1f}) on f = x^2 + 3y^2, where f = {S.bowl(p):.6f}")
print(f"    a step of {step} ALONG the contour direction changes f by "
      f"{S.bowl(p + step * unit(along)) - S.bowl(p):+.3e}")
print(f"    the same step ACROSS it, up the gradient, changes f by "
      f"{S.bowl(p + step * unit(g)) - S.bowl(p):+.3e}")
print(f"    the ratio of those two changes is about "
      f"{abs((S.bowl(p + step * unit(g)) - S.bowl(p)) / (S.bowl(p + step * unit(along)) - S.bowl(p))):.0f} to 1")
gain_across = S.bowl(p + step * unit(g)) - S.bowl(p)
gain_along = S.bowl(p + step * unit(along)) - S.bowl(p)
assert abs(gain_along) < 1e-5
assert gain_across > 100 * abs(gain_along)
print()
print(f"    and the gradient's length says how fast: {magnitude(g):.6f} units of f")
print(f"    per unit of distance, so a step of {step} up the gradient should gain")
print(f"    about {magnitude(g) * step:.6f}, against the {gain_across:.6f} actually measured.")
assert abs(gain_across - magnitude(g) * step) < 1e-4

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

"""Constant gradients, bowl gradients, and the three faces of a zero gradient.

Run from inside `examples/`:

    ../.venv/bin/python3 05_flat_ground_three_ways.py
"""

from __future__ import annotations

import numpy as np

import surfaces as S
from gradients import angle_degrees, gradient, magnitude, unit

print(__doc__.splitlines()[0])
print()

# --------------------------------------------------------------------------
print("1. A plane: the same gradient everywhere, however far you walk")
# --------------------------------------------------------------------------
#
# f(x, y) = 3x - 2y + 5 is a flat tilted sheet. Its partial derivatives are 3
# and -2, and neither one mentions x or y, so the gradient is the constant
# vector [3, -2]. A tilted sheet has one slope and one uphill direction and
# they are the same at every point on it -- which is why a linear model has
# nothing to optimise: there is no bottom to fall into.

print("  f(x, y) = 3x - 2y + 5,  grad f = (3, -2) everywhere")
print()
print(f"  {'point':>18}  {'f':>12}  {'gradient':>22}  {'|gradient|':>12}  {'bearing':>10}")
seen = []
for p in ((0.0, 0.0), (1.0, 1.0), (-40.0, 17.5), (0.001, 0.002)):
    g = gradient(S.plane, p)
    seen.append(g)
    shown = "[" + ", ".join(f"{v:9.6f}" for v in g) + "]"
    print(f"  {str(p):>18}  {S.plane(p):12.4f}  {shown:>22}  {magnitude(g):12.7f}"
          f"  {angle_degrees(g):9.3f}d")
    assert abs(g[0] - 3.0) < S.GRADIENT_TOL
    assert abs(g[1] + 2.0) < S.GRADIENT_TOL

spread = float(np.max(np.abs(np.array(seen) - np.array(seen[0]))))
print()
print(f"  Largest disagreement between any two of those gradients: {spread:.3e}")
print("  -- which is floating-point noise, not variation. The value of f")
print("  swings from +6 to -150 across those points; the gradient does not")
print("  change at all.")
assert spread < S.GRADIENT_TOL

# --------------------------------------------------------------------------
print()
print("1b. Where that stops being measurable, and exactly why")
# --------------------------------------------------------------------------
#
# This section was not planned. It was found by putting (1000, -1000) in the
# table above and watching the assertion fail, and it is kept because it is
# more useful than the tidy version would have been.
#
# The gradient of a plane is constant, so a numerical estimate ought to be
# equally good anywhere. It is not. The central difference computes
#
#     ( f(x+h) - f(x-h) ) / 2h
#
# and the two values of f are stored with a relative error of about one
# machine epsilon EACH. At (1000, -1000) the function is worth about 5005, so
# each stored value carries an absolute error of roughly 5005 * 2.2e-16, and
# the subtraction keeps that error while the division by 2h = 2e-5 multiplies
# it by fifty thousand. Predicted noise: eps * |f| / (2h).
#
# That is not an approximation of the trouble; it is the trouble, and the
# prediction can be checked.

eps = float(np.finfo(float).eps)
print("  The same constant gradient, measured further and further from the origin:")
print()
print(f"  {'point':>26}  {'|f|':>13}  {'measured df/dx':>16}  {'error':>11}  {'eps|f|/2h':>11}")
for p in ((1.0, 1.0), (-40.0, 17.5), (1000.0, -1000.0),
          (100000.0, -100000.0), (10000000.0, -10000000.0)):
    g = gradient(S.plane, p)
    error = float(np.max(np.abs(g - np.array([3.0, -2.0]))))
    predicted = eps * abs(S.plane(p)) / (2.0 * S.H_DEFAULT)
    print(f"  {str(p):>26}  {abs(S.plane(p)):13.1f}  {g[0]:16.10f}  {error:11.3e}"
          f"  {predicted:11.3e}")
    # The prediction is a bound with a small constant, not an identity: the
    # measured error must be of the same order, never wildly above it.
    assert error < 3.0 * predicted + 1e-12, (p, error, predicted)

print()
print("  The last two columns track each other across seven orders of")
print("  magnitude. By ten million the estimate of a gradient that is exactly")
print("  3 has lost its fourth decimal place, and nothing about the calculus")
print("  went wrong -- only the arithmetic.")
print()
print("  This is the boundary on everything else in the lab. The tolerance of")
print(f"  {S.GRADIENT_TOL:g} that every other assertion uses is only achievable because every")
print("  probe point in `surfaces.py` keeps |f| small. Feed a numerical")
print("  gradient a loss of a hundred thousand and it will hand you a")
print("  confident answer with four good digits in it. Autodiff, which")
print("  differentiates the expression rather than sampling it, does not have")
print("  this failure mode at all -- which is the first of the two reasons")
print("  nobody trains a model this way.")

# --------------------------------------------------------------------------
print()
print("2. A bowl: the gradient points AWAY from the minimum, and grows with distance")
# --------------------------------------------------------------------------
#
# f(x, y) = x^2 + 3y^2 has its minimum at the origin. The gradient points
# uphill, and uphill from anywhere on a bowl is away from the bottom. Two
# consequences that Day 111 depends on: the NEGATIVE gradient always points
# roughly back towards the minimum, and it gets shorter as you approach, so
# the steps naturally shrink near the answer.

print(f"  {'point':>16}  {'distance from origin':>21}  {'gradient':>22}"
      f"  {'|gradient|':>12}  {'points away?':>13}")
for p in ((0.5, 0.5), (1.0, 1.0), (2.0, 2.0), (4.0, 4.0), (-3.0, 1.0)):
    q = np.array(p)
    g = gradient(S.bowl, p)
    outward = float(np.dot(unit(g), unit(q)))   # positive means "away from origin"
    shown = "[" + ", ".join(f"{v:9.5f}" for v in g) + "]"
    print(f"  {str(p):>16}  {magnitude(q):21.6f}  {shown:>22}  {magnitude(g):12.6f}"
          f"  {'yes' if outward > 0 else 'no':>13}")
    assert outward > 0.0

print()
print("  Every row points away from the bottom, and the length grows with")
print("  distance. Walk against it and you head back down -- which is the")
print("  entire algorithm of Day 111, and the reason the steps get smaller by")
print("  themselves as the answer gets closer.")

print()
print("  It is not, however, aimed exactly at the origin, because the bowl is")
print("  elliptical. Compare the gradient's bearing with the bearing straight")
print("  back to the minimum:")
print(f"    {'point':>16}  {'bearing of -gradient':>21}  {'bearing to origin':>19}  {'off by':>9}")
for p in ((1.0, 1.0), (3.0, 0.5), (0.5, 3.0)):
    g = gradient(S.bowl, p)
    back = angle_degrees(-g)
    straight = angle_degrees(-np.array(p))
    print(f"    {str(p):>16}  {back:20.3f}d  {straight:18.3f}d  {abs(back - straight):8.3f}d")

print()
print("  On a circular bowl those two would agree exactly. On this one they do")
print("  not, and that mismatch is precisely what makes gradient descent")
print("  zig-zag down a narrow valley instead of walking straight in.")

# --------------------------------------------------------------------------
print()
print("3. Zero gradient, three completely different points")
# --------------------------------------------------------------------------
#
# All three surfaces have gradient [0, 0] at the origin. The gradient is
# identical. The points are not. This is the honest limit of everything the
# day has built: a zero gradient tells you the ground is level, and stops.

print("  Walk 0.1 in eight directions from the origin and record what f does:")
print()
radius = 0.1
bearings = np.arange(0, 360, 45)
print(f"    {'surface':>9}  " + "  ".join(f"{b:>7}d" for b in bearings) + "   verdict")
for name, kind, why in S.STATIONARY_AT_ORIGIN:
    f = S.SURFACES[name][0]
    g = gradient(f, (0.0, 0.0))
    assert magnitude(g) < S.GRADIENT_TOL
    changes = []
    for b in bearings:
        a = np.radians(float(b))
        step = radius * np.array([np.cos(a), np.sin(a)])
        changes.append(f(step) - f(np.array([0.0, 0.0])))
    row = "  ".join(f"{c:+8.4f}" for c in changes)
    up = sum(1 for c in changes if c > 1e-12)
    down = sum(1 for c in changes if c < -1e-12)
    if down == 0:
        verdict = "all up -- a minimum"
    elif up == 0:
        verdict = "all down -- a maximum"
    else:
        level = len(changes) - up - down
        verdict = f"{up} up, {down} down, {level} flat -- a saddle"
    print(f"    {name:>9}  {row}   {verdict}")
    assert kind in verdict

print()
print("  The saddle's four flat entries are not rounding: on x^2 - y^2 the")
print("  diagonals are exactly where x^2 equals y^2, so f is unchanged along")
print("  them. They are the two contour lines that cross AT the saddle, which")
print("  is what makes a saddle a saddle.")
print()
print("  Same gradient. Three different answers. Nothing in the gradient")
print("  distinguishes them, and no amount of care computing it would help,")
print("  because the information simply is not in there: the gradient is built")
print("  from FIRST derivatives, and which kind of stationary point this is")
print("  depends on the SECOND ones. That object has a name -- the Hessian --")
print("  and this course does not develop it here.")

# --------------------------------------------------------------------------
print()
print("4. Why the saddle is the one that matters")
# --------------------------------------------------------------------------
#
# In two dimensions a saddle is a curiosity. In the parameter space of a
# model, where there are millions of directions rather than two, a point where
# every partial derivative is zero is overwhelmingly more likely to be a
# saddle than a true minimum -- because being a minimum requires the surface
# to curve upward in EVERY one of those millions of directions at once, and
# being a saddle only requires one direction to disagree.

print("  On f(x, y) = x^2 - y^2 the origin is level, and the two axes disagree:")
for label, direction in (("along x", (1.0, 0.0)), ("along y", (0.0, 1.0))):
    for r in (0.1, 0.5, 1.0):
        step = r * np.array(direction)
        print(f"    {label}, distance {r:>4}:  f = {S.saddle(step):+8.4f}")
    print()

print("  Walking east from the origin, f rises. Walking north, it falls. The")
print("  gradient at the origin is zero in both cases, and an optimiser that")
print("  stops when the gradient is zero would stop right there, in a place it")
print("  could have escaped by moving a millimetre north.")
assert S.saddle((0.5, 0.0)) > 0
assert S.saddle((0.0, 0.5)) < 0
assert magnitude(gradient(S.saddle, (0.0, 0.0))) < S.GRADIENT_TOL

print()
print("  And the gradient near the saddle is small without being zero, which")
print("  is the practical problem: progress crawls rather than stopping,")
print("  and it is hard to tell the two apart from the outside.")
print(f"    {'distance from the saddle':>26}  {'|gradient|':>12}")
for r in (1.0, 0.1, 0.01, 0.001):
    p = (r, r)
    print(f"    {r:26}  {magnitude(gradient(S.saddle, p)):12.6f}")

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

"""Choosing h, Day 108's U-curve in two dimensions, and what numpy.gradient does.

Run from inside `examples/`:

    ../.venv/bin/python3 06_step_size_and_the_u_curve.py
"""

from __future__ import annotations

import numpy as np

import surfaces as S
from gradients import forward_partial, gradient, partial

print(__doc__.splitlines()[0])
print()

POINT = (2.0, 1.0)
EXACT_DX = float(S.cubic_gradient(POINT)[0])

# --------------------------------------------------------------------------
print("1. On a cubic, the truncation error is not merely small -- it is exactly h^2")
# --------------------------------------------------------------------------
#
# Script 01 showed that a central difference is algebraically exact on a
# quadratic. f(x, y) = x^3 + x*y^2 is the first surface here that is not, and
# its error can be written down in closed form rather than bounded:
#
#     ((x+h)^3 - (x-h)^3) / (2h)
#   = (x^3 + 3x^2h + 3xh^2 + h^3 - x^3 + 3x^2h - 3xh^2 + h^3) / (2h)
#   = (6x^2 h + 2h^3) / (2h)
#   = 3x^2 + h^2
#
# and the exact partial is 3x^2. So the numerical answer overshoots by exactly
# h squared, with no other terms at all. The y^2 x term contributes nothing to
# the x-partial's error because it is linear in x.

print(f"  f(x, y) = x^3 + x*y^2 at {POINT}. Exact df/dx = 3x^2 + y^2 = {EXACT_DX}")
print()
print(f"  {'h':>10}  {'numerical df/dx':>20}  {'error':>16}  {'h^2':>16}  {'relative gap':>13}")
for k in (1, 2, 3):
    h = 10.0 ** (-k)
    value = partial(S.cubic, POINT, 0, h)
    error = value - EXACT_DX
    gap = abs(error - h * h) / (h * h)
    print(f"  {h:10.0e}  {value:20.14f}  {error:16.12f}  {h * h:16.12f}  {gap:13.3e}")
    assert gap < 1e-5, (h, gap)

print()
print("  The error column and the h^2 column are the same column. This is the")
print("  clearest statement available of what 'second-order accurate' means:")
print("  divide the step by ten and the method error divides by a hundred.")

# --------------------------------------------------------------------------
print()
print("2. Day 108's U-curve, on a partial derivative")
# --------------------------------------------------------------------------
#
# Shrinking h forever does not work, and Day 108 already showed why in one
# dimension: the method error falls but the ROUNDOFF error rises, because
# subtracting two numbers that are nearly equal throws away the leading digits
# they had in common. Nothing about that changes when the function has more
# than one input.
#
# Two curves are printed, because the choice between them is the whole reason
# this lab uses a central difference.

print(f"  {'h':>10}  {'central error':>16}  {'forward error':>16}   shape")
central = {}
forward = {}
for k in range(0, 15):
    h = 10.0 ** (-k)
    c = abs(partial(S.cubic, POINT, 0, h) - EXACT_DX)
    f = abs(forward_partial(S.cubic, POINT, 0, h) - EXACT_DX)
    central[h] = c
    forward[h] = f
    bar = "#" * max(0, int(round(16 + np.log10(max(c, 1e-16)))))
    print(f"  {h:10.0e}  {c:16.3e}  {f:16.3e}   {bar}")

best_central = min(central, key=central.get)
best_forward = min(forward, key=forward.get)
eps = float(np.finfo(float).eps)
print()
print(f"  best h for the central difference: {best_central:.0e}"
      f"   (error {central[best_central]:.3e})")
print(f"  best h for the forward difference: {best_forward:.0e}"
      f"   (error {forward[best_forward]:.3e})")
print()
print("  Theory says the trough sits where the two error sources balance:")
print(f"    central: around the cube root of machine epsilon = {eps ** (1 / 3):.3e}")
print(f"    forward: around the square root of machine epsilon = {eps ** 0.5:.3e}")
print()
print("  Both predictions land on the measured trough to within one decade.")
assert best_central == 1e-05
assert best_forward == 1e-08
assert central[best_central] < forward[best_forward]

print()
print("  Read the two error columns at h = 1e-5, the step this lab uses:")
print(f"    central  {central[1e-05]:.3e}")
print(f"    forward  {forward[1e-05]:.3e}")
print(f"    the central difference is {forward[1e-05] / central[1e-05]:,.0f} times more accurate here,")
print("    for one extra evaluation of f per input. That is the trade, and it")
print("    is not close.")
print()
print("  And the far end of the table is the part worth remembering: at")
print(f"  h = 1e-14 the central difference is out by {central[1e-14]:.3e}, which is")
print(f"  {central[1e-14] / central[1e-01]:.0f} times WORSE than the answer at h = 0.1 -- a step a trillion")
print("  times bigger. Shrinking h past the trough does not buy a slightly")
print("  worse answer. It buys nonsense, confidently.")
assert central[1e-14] > central[1e-01]
assert central[1e-05] < central[1e-01]

# --------------------------------------------------------------------------
print()
print("3. What happens to the whole gradient, not just one partial")
# --------------------------------------------------------------------------

print(f"  {'h':>10}  {'gradient of the cubic at (2, 1)':>36}  {'max error':>12}")
for k in (1, 3, 5, 8, 12):
    h = 10.0 ** (-k)
    g = gradient(S.cubic, POINT, h)
    err = float(np.max(np.abs(g - S.cubic_gradient(POINT))))
    shown = "[" + ", ".join(f"{v:16.10f}" for v in g) + "]"
    print(f"  {h:10.0e}  {shown:>36}  {err:12.3e}")

print()
print("  Both components degrade together, because both are computed the same")
print("  way. There is no step size that is right for one and wrong for the")
print("  other here -- though on a function whose inputs have wildly different")
print("  scales there would be, which is an argument for scaling your inputs")
print("  before you differentiate anything.")

# --------------------------------------------------------------------------
print()
print("4. numpy.gradient does something related and different")
# --------------------------------------------------------------------------
#
# NumPy has a function called `gradient`, and reaching for it here would be a
# mistake -- not because it is bad, but because it answers a different
# question. Ours takes a FUNCTION and a point. NumPy's takes an ARRAY of
# values already sampled on a grid, and returns the differences between
# neighbouring samples. It cannot be asked for the gradient at a point that is
# not a grid point, and it cannot choose its own step, because the step is
# whatever spacing the data already has.

xs = np.linspace(0.0, 4.0, 9)
ys = np.linspace(0.0, 4.0, 9)
X, Y = np.meshgrid(xs, ys, indexing="ij")
spacing = float(xs[1] - xs[0])

print(f"  A 9 by 9 grid over [0, 4] x [0, 4], spacing {spacing}.")
print()
print("  On the bowl x^2 + 3y^2, whose gradient a central difference gets")
print("  exactly right at any step, numpy.gradient is exact in the interior:")
Z = X * X + 3.0 * Y * Y
gx, gy = np.gradient(Z, xs, ys)
print(f"    interior sample at (x, y) = ({xs[2]}, {ys[2]}):"
      f" numpy [{gx[2, 2]:.6f}, {gy[2, 2]:.6f}]"
      f"  exact [{S.bowl_gradient((xs[2], ys[2]))[0]:.6f},"
      f" {S.bowl_gradient((xs[2], ys[2]))[1]:.6f}]")
assert abs(gx[2, 2] - 2.0 * xs[2]) < 1e-12
assert abs(gy[2, 2] - 6.0 * ys[2]) < 1e-12

print()
print("  but not at the edge, because by default it drops to a one-sided")
print("  first-order formula there:")
print(f"    corner sample at (0.0, 0.0):  numpy [{gx[0, 0]:.6f}, {gy[0, 0]:.6f}]"
      f"   exact [0.0, 0.0]")
assert abs(gx[0, 0] - 0.5) < 1e-12
assert abs(gy[0, 0] - 1.5) < 1e-12

gx2, gy2 = np.gradient(Z, xs, ys, edge_order=2)
print(f"    the same corner with edge_order=2: [{gx2[0, 0]:.6f}, {gy2[0, 0]:.6f}]"
      "   exact, this time")
assert abs(gx2[0, 0]) < 1e-12 and abs(gy2[0, 0]) < 1e-12
print()
print("  That default is worth knowing about before it costs you an afternoon:")
print("  every interior value is second-order accurate and every boundary")
print("  value is first-order, unless you ask otherwise.")

print()
print("  On the cubic, where the method error is real, the difference between")
print("  the two functions becomes the point:")
Zc = X ** 3 + X * Y * Y
cgx, _cgy = np.gradient(Zc, xs, ys, edge_order=2)
i = j = 4
p = (float(xs[i]), float(ys[j]))
exact = float(S.cubic_gradient(p)[0])
ours = float(partial(S.cubic, p, 0))
print(f"    at (x, y) = {p}, exact df/dx = {exact}")
print(f"      numpy.gradient on the sampled array : {cgx[i, j]:14.10f}"
      f"   error {abs(cgx[i, j] - exact):.3e}")
print(f"      grid spacing squared                 : {spacing ** 2:14.10f}")
print(f"      our gradient, on the function itself : {ours:14.10f}"
      f"   error {abs(ours - exact):.3e}")
assert abs(abs(cgx[i, j] - exact) - spacing ** 2) < 1e-12
assert abs(ours - exact) < S.GRADIENT_TOL

print()
print("  It is the same h^2 law from section 1 -- but h is now the grid")
print("  spacing, which is fixed by the data you were given. You cannot")
print("  shrink it without going back and sampling more finely, and if the")
print("  samples came from a sensor you may not be able to at all.")
print()
print("  So: numpy.gradient when you HAVE an array of values -- an image, a")
print("  height field, a measured series. Our `gradient` when you have a")
print("  function you can call. They are not competitors.")

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

"""Every parameter of a model gets a partial derivative. That collection is the gradient.

Run from inside `examples/`:

    ../.venv/bin/python3 07_one_partial_per_parameter.py
"""

from __future__ import annotations

import numpy as np

import surfaces as S
from gradients import gradient, magnitude, unit

print(__doc__.splitlines()[0])
print()

# --------------------------------------------------------------------------
print("1. The smallest thing that is honestly a model")
# --------------------------------------------------------------------------
#
# Three parameters and four invented samples. `pred = w1*a + w2*b + c`, and
# the loss is the mean of the squared differences between what the model says
# and what the target says.
#
# Nothing about this is a toy version of the idea. It IS the idea. A network
# with a hundred million parameters differs from this in exactly one respect:
# the number of parameters.

print("  Four invented samples. Nothing here was measured; the numbers were")
print("  chosen so every step below can be checked with a pencil.")
print()
print(f"    {'a':>6}  {'b':>6}  {'target':>8}")
for a, b, target in S.SAMPLES:
    print(f"    {a:6.1f}  {b:6.1f}  {target:8.1f}")

params = S.START_PARAMS
w1, w2, c = params
print()
print(f"  Parameters to start with: w1 = {w1}, w2 = {w2}, c = {c}")
print()
print(f"    {'a':>6}  {'b':>6}  {'prediction':>11}  {'target':>8}  {'residual':>9}  {'squared':>9}")
total = 0.0
for a, b, target in S.SAMPLES:
    pred = w1 * a + w2 * b + c
    residual = pred - target
    total += residual * residual
    print(f"    {a:6.1f}  {b:6.1f}  {pred:11.2f}  {target:8.1f}  {residual:9.2f}  {residual ** 2:9.2f}")
print(f"    {'':>6}  {'':>6}  {'':>11}  {'':>8}  {'sum':>9}  {total:9.2f}")
print(f"    mean squared error over {len(S.SAMPLES)} samples: {total / len(S.SAMPLES)}")
assert S.model_loss(params) == total / len(S.SAMPLES)
assert S.model_loss(params) == 22.5

# --------------------------------------------------------------------------
print()
print("2. One partial derivative per parameter")
# --------------------------------------------------------------------------
#
# The loss is a function of THREE inputs -- and they are not the data. The
# data is fixed; what varies, and what the derivative is taken with respect
# to, is the parameters. That swap is the thing worth stopping on: a and b are
# constants inside this function, and w1, w2 and c are the variables.

numeric = gradient(S.model_loss, params)
exact = S.model_loss_gradient(params)
print("  Nudge each parameter on its own, hold the other two still:")
print()
print(f"    {'parameter':>10}  {'numerical':>18}  {'exact, by hand':>15}  {'error':>11}")
for i, label in enumerate(("w1", "w2", "c")):
    print(f"    {label:>10}  {numeric[i]:18.12f}  {exact[i]:15.4f}"
          f"  {abs(numeric[i] - exact[i]):11.3e}")
    assert abs(numeric[i] - exact[i]) < S.GRADIENT_TOL

print()
print("  Working the first one by hand, to show there is no magic:")
print("    L = (1/4) sum (w1*a + w2*b + c - y)^2")
print("    dL/dw1 = (2/4) sum (w1*a + w2*b + c - y) * a       [chain rule, Day 110]")
print("           = 0.5 * ( -4*1  +  -3*2  +  -8*3  +  -1*0 )")
print("           = 0.5 * (-34)")
print("           = -17")
assert exact[0] == -17.0
assert exact[1] == -18.0
assert exact[2] == -8.0

print()
print(f"  So grad L = [{exact[0]:.0f}, {exact[1]:.0f}, {exact[2]:.0f}] -- three numbers, one per parameter,")
print("  and they are all negative, which says every parameter is currently")
print("  too small: increasing any of them increases the loss's rate of")
print("  DEcrease. Day 111 acts on that.")

# --------------------------------------------------------------------------
print()
print("3. The gradient is still just a vector, so everything from Day 99 applies")
# --------------------------------------------------------------------------

g = exact
print(f"  gradient        [{g[0]:.0f}, {g[1]:.0f}, {g[2]:.0f}]")
print(f"  length          {magnitude(g):.6f}   -- how steep the loss surface is here")
u = unit(g)
print(f"  unit gradient   [{u[0]:.6f}, {u[1]:.6f}, {u[2]:.6f}]")
print(f"  its length      {magnitude(u):.15f}")
assert abs(magnitude(u) - 1.0) < 1e-12
print()
print("  There is no picture of this one. The input space has three dimensions")
print("  and the surface would need four, which nobody can draw. Every single")
print("  statement from the two-dimensional case survives the move anyway:")
print("  the gradient points the steepest way up, its length says how steep,")
print("  and it is perpendicular to the level set -- which is now a surface")
print("  rather than a curve. That transfer is the reason the day was spent on")
print("  pictures of hills.")

# --------------------------------------------------------------------------
print()
print("4. One step against the gradient, to show it is not a claim")
# --------------------------------------------------------------------------
#
# Day 111 is the day this becomes an algorithm. But the single step can be
# taken here in three lines, and it is more convincing than any amount of
# assurance that it would work.

print(f"  {'step size':>10}  {'new parameters':>34}  {'loss':>14}  {'change':>12}")
before = S.model_loss(params)
print(f"  {'0 (start)':>10}  "
      f"{'[' + ', '.join(f'{v:9.5f}' for v in params) + ']':>34}  {before:14.8f}  {'':>12}")
improved = 0
steps = (0.001, 0.005, 0.01, 0.02, 0.05, 0.1, 0.15, 0.2)
losses = {}
for step in steps:
    moved = np.array(params) - step * np.array(exact)
    after = S.model_loss(moved)
    shown = "[" + ", ".join(f"{v:9.5f}" for v in moved) + "]"
    print(f"  {step:10.3f}  {shown:>34}  {after:14.8f}  {after - before:+12.6f}")
    losses[step] = after
    if after < before:
        improved += 1

print()
best_step = min(losses, key=losses.get)
print(f"  {improved} of the {len(steps)} step sizes reduced the loss, and the best of them was")
print(f"  {best_step}, which brought it from {before} down to {losses[best_step]:.5f}.")
print(f"  Past that the loss climbs again -- {losses[0.2]:.5f} at a step of 0.2, which is")
print("  WORSE than where it started. So stepping against the gradient is the")
print("  right DIRECTION, and how far to go along it is a separate question")
print("  with its own name -- the learning rate -- and its own way of going")
print("  wrong. Day 111 is about both.")
assert improved >= 4
assert losses[0.2] > before
assert losses[best_step] < before
assert S.model_loss(np.array(params) - 0.001 * np.array(exact)) < before

# --------------------------------------------------------------------------
print()
print("5. The cost, and why nobody trains a model this way")
# --------------------------------------------------------------------------
#
# `gradient` calls `partial` once per input, and `partial` evaluates f twice.
# So a numerical gradient of a function of n inputs costs 2n evaluations of
# the whole function. That is the entire argument for automatic
# differentiation, and it is arithmetic rather than opinion.

class Counter:
    """Wraps a function and counts how many times it is actually called."""

    def __init__(self, f):
        self.f = f
        self.calls = 0

    def __call__(self, point):
        self.calls += 1
        return self.f(point)


counted = Counter(S.model_loss)
gradient(counted, params)
print(f"  Parameters: {len(params)}")
print(f"  Evaluations of the loss to get one gradient: {counted.calls}")
print(f"  Which is 2 per parameter: {2 * len(params)}")
assert counted.calls == 2 * len(params)

print()
print("  Now scale that. One forward pass of a model is one evaluation of the")
print("  loss, so the table below is in units of 'complete forward passes")
print("  through the entire network, over the entire batch, per single")
print("  training step':")
print()
print(f"    {'parameters':>14}  {'forward passes for ONE numerical gradient':>44}")
for n in (3, 1_000, 1_000_000, 1_000_000_000):
    print(f"    {n:>14,}  {2 * n:>44,}")

print()
print("  A million-parameter model would need two million forward passes to")
print("  take one step. Reverse-mode automatic differentiation -- what")
print("  PyTorch's autograd and JAX's grad do -- gets the whole gradient for a")
print("  cost of roughly ONE forward pass plus one backward pass, no matter")
print("  how many parameters there are, and gets it exactly rather than to")
print("  within h^2. That is not an optimisation of the method in this file.")
print("  It is a different method, and it is the reason training large models")
print("  is possible at all.")
print()
print("  None of those libraries is installed in this lab and no output from")
print("  them is reproduced anywhere in it. What numerical differentiation IS")
print("  still good for is checking one: if your hand-written backward pass")
print("  disagrees with a numerical gradient on a small example, the")
print("  hand-written one is wrong. That check has a name -- gradient")
print("  checking -- and this file is a working implementation of it.")

print()
print("07_one_partial_per_parameter.py: every assertion held.")

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

Both `examples/` and `starter/` contain modules called `gradients` and
`surfaces`, 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 `gradients` 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
`gradients`, `surfaces` 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 ("gradients", "surfaces", "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]

"""Numerical partial derivatives and gradients, built from nothing.

Eleven functions. The first two are the whole day: `partial` asks what happens
when you nudge ONE input and hold the rest still, and `gradient` collects one
of those per input into a vector. Everything else in this file is built out of
those two.

Nothing here knows anything about the functions it is handed. It cannot see
inside them, cannot differentiate them symbolically, and does not try. It
evaluates them at points and subtracts. That is both the strength (it works on
anything you can call) and the weakness (it costs two evaluations per input,
which is why nobody trains a neural network this way -- see the lesson's
alternatives section, and Day 110).
"""

from __future__ import annotations

import numpy as np

from surfaces import H_DEFAULT


def partial(f, point, index, h=H_DEFAULT):
    """The partial derivative of `f` with respect to input number `index`.

    Nudge coordinate `index` up by h and down by h, leave every other
    coordinate exactly where it was, and divide the change in f by the total
    distance moved, 2h. That is Day 108's central difference with one word
    added: "and hold everything else still".

        df/dx_i  ~  ( f(... x_i + h ...) - f(... x_i - h ...) ) / (2h)

    The forward difference ( f(x+h) - f(x) ) / h would also work and costs one
    fewer evaluation when you already have f(x). It is also markedly worse:
    its error shrinks like h where the central difference's shrinks like h
    squared. Script 06 measures both.
    """
    base = np.asarray(point, dtype=float)
    up = base.copy()
    down = base.copy()
    up[index] += h
    down[index] -= h
    return float((f(up) - f(down)) / (2.0 * h))


def forward_partial(f, point, index, h=H_DEFAULT):
    """The same thing with a one-sided step, kept only for the comparison."""
    base = np.asarray(point, dtype=float)
    up = base.copy()
    up[index] += h
    return float((f(up) - f(base)) / h)


def gradient(f, point, h=H_DEFAULT):
    """The gradient: one partial derivative per input, collected into a vector.

    This is the definition and there is nothing hidden in it. If the point has
    two coordinates you get a vector of two numbers; if it has three you get
    three; if it has a million you get a million, and the loop below is why
    nobody does it that way.

    Written with the nabla symbol as grad f, and it is a VECTOR living in the
    same space as the input, not a number and not a point on the surface.
    """
    base = np.asarray(point, dtype=float)
    return np.array([partial(f, base, i, h) for i in range(base.size)])


def magnitude(vector):
    """The length of a vector -- Day 99's Euclidean norm, unchanged.

    Applied to a gradient it answers "how steep is the steepest way up",
    in units of f per unit of distance travelled in the input space.
    """
    return float(np.sqrt(np.dot(np.asarray(vector, dtype=float),
                                np.asarray(vector, dtype=float))))


def unit(vector):
    """The same direction, scaled to length 1.

    A direction has to be a unit vector before a directional derivative means
    anything: without that, doubling the vector would double the answer and
    the "rate of change in this direction" would depend on how long an arrow
    you happened to draw.
    """
    v = np.asarray(vector, dtype=float)
    length = magnitude(v)
    if length == 0.0:
        raise ValueError("the zero vector has no direction")
    return v / length


def directional_derivative(f, point, direction, h=H_DEFAULT):
    """How fast f changes if you walk from `point` along `direction`.

    Two lines, and the second is Day 103 doing real work: normalise the
    direction, then dot it with the gradient. That the dot product is the
    right answer is not obvious and is not asserted here -- script 03 checks
    it against a direct measurement, which is a straight central difference
    taken ALONG the direction rather than along an axis.
    """
    u = unit(direction)
    return float(np.dot(gradient(f, point, h), u))


def directional_derivative_direct(f, point, direction, h=H_DEFAULT):
    """The same quantity measured without ever forming a gradient.

    Step h forward along the direction and h back along it, and divide by 2h.
    No partials, no dot product, no assumption that the two agree. This
    function exists to check the one above.
    """
    base = np.asarray(point, dtype=float)
    u = unit(direction)
    return float((f(base + h * u) - f(base - h * u)) / (2.0 * h))


def sweep_directions(f, point, n=None, h=H_DEFAULT):
    """Try n directions evenly spaced around the circle; report each rate.

    Returns (angles_in_radians, rates). Only meaningful for a function of two
    inputs, which is the case the reader can draw.
    """
    from surfaces import N_DIRECTIONS

    if n is None:
        n = N_DIRECTIONS
    angles = np.linspace(0.0, 2.0 * np.pi, n, endpoint=False)
    rates = np.array([
        directional_derivative_direct(f, point, np.array([np.cos(a), np.sin(a)]), h)
        for a in angles
    ])
    return angles, rates


def angle_degrees(vector):
    """The compass bearing of a 2-D vector, measured from the positive x-axis,
    reported in [0, 360) degrees so two angles can be compared directly."""
    v = np.asarray(vector, dtype=float)
    return float(np.degrees(np.arctan2(v[1], v[0])) % 360.0)


def angular_gap_degrees(a, b):
    """The smaller of the two ways round between two bearings, in degrees.

    Without the wrap-around, 359.6 and 0.1 would look 359.5 degrees apart
    instead of 0.5, and the steepest-ascent check would fail for a reason
    that has nothing to do with calculus.
    """
    raw = abs(a - b) % 360.0
    return float(min(raw, 360.0 - raw))


def contour_chord(f, contour, level, t, delta):
    """A unit vector along the contour of f, built without using the gradient.

    Takes two points on the exact algebraic contour, at parameters t and
    t + delta, and returns the unit vector from the first to the second, plus
    the two points. As delta shrinks the chord approaches the tangent, and the
    tangent is the thing the gradient is claimed to be perpendicular to.

    The value of f at both points is returned as well, so the caller can check
    that the parametrisation really does stay on one level rather than trust
    the algebra.
    """
    p = contour(level, t)
    q = contour(level, t + delta)
    return unit(q - p), p, q, f(p), f(q)

"""The surfaces this lab measures, and the exact gradients to check against.

Every function here is a function of several inputs. Every one has a gradient
you can work out with a pencil in under a minute, which is the entire point:
the numerical machinery in `gradients.py` is only trustworthy if there is
something exact to hold it against.

All the data is invented. None of it is a measurement of anything real. What
IS real is every number the lab prints about these functions, because those are
computed from these definitions at run time.

Read this file. Do not change it -- the tests compare against the values
written down here.
"""

from __future__ import annotations

import numpy as np

# --------------------------------------------------------------------------
# The step size, and why this one
# --------------------------------------------------------------------------
#
# Day 108 established the shape: a central difference has a truncation error
# that shrinks like h squared and a roundoff error that GROWS like 1/h, so the
# total error is U-shaped in h and the best step is somewhere in the middle.
# For a central difference on float64 the trough sits near the cube root of
# the machine epsilon, which is about 6e-6. This lab uses 1e-5, which is close
# enough to the bottom of that trough to be within a factor of two of the best
# achievable error on every surface here, and is a round number a reader can
# remember. Script 06 sweeps h over twelve orders of magnitude and prints the
# curve rather than asking you to take this on trust.
H_DEFAULT = 1.0e-5

# The tolerance every gradient assertion in this lab uses.
#
# It is set from what the arithmetic can achieve, not from what makes the
# tests pass. Four of the six surfaces below are at most quadratic in each
# variable, and a central difference is ALGEBRAICALLY EXACT for those -- the
# h-squared term is multiplied by a third derivative that is zero, so the only
# error left is floating-point roundoff, which lands around 1e-11 at h = 1e-5
# for points of this size. The one genuinely cubic surface has a truncation
# error of exactly h squared, which is 1e-10 here. GRADIENT_TOL is three
# orders of magnitude above the worst of those, which leaves room for a
# different processor's rounding without leaving room for a wrong answer.
GRADIENT_TOL = 1.0e-8

# Directions sampled around the full circle when the lab checks that the
# gradient really is the steepest way up. With 360 evenly spaced directions
# the closest sample sits within 0.5 degrees of any given angle, so the
# tolerance below is that bound with a little slack for the wrap-around at
# 360 degrees.
N_DIRECTIONS = 360
ANGLE_TOL_DEGREES = 1.0

# The step taken along a contour when the lab checks perpendicularity. The
# chord between two points on a curve is not the tangent; it differs from it
# by an angle of order delta, so the dot product of a unit gradient with a
# unit chord is of order delta rather than exactly zero. Script 04 halves
# delta four times and prints the dot product each time, so the reader watches
# it shrink instead of being handed a tolerance.
CONTOUR_DELTA = 1.0e-5
CONTOUR_DOT_TOL = 1.0e-4

# Used where a random point or a random direction is wanted. Seeded, so every
# number in `expected-output/` is reproducible.
SEED = 109


# --------------------------------------------------------------------------
# 1. A quadratic bowl -- the shape every optimisation picture is drawn on
# --------------------------------------------------------------------------

def bowl(point):
    """f(x, y) = x^2 + 3y^2. A bowl with its lowest point at the origin.

    The 3 makes it an elliptical bowl rather than a circular one: it is three
    times steeper in y than in x at the same distance out, which is exactly
    the situation that makes gradient descent zig-zag on Day 111.
    """
    x, y = point
    return x * x + 3.0 * y * y


def bowl_gradient(point):
    """grad f = (2x, 6y). Differentiate x^2 + 3y^2 one variable at a time."""
    x, y = point
    return np.array([2.0 * x, 6.0 * y])


# --------------------------------------------------------------------------
# 2. A plane -- the function whose gradient is the same everywhere
# --------------------------------------------------------------------------

def plane(point):
    """f(x, y) = 3x - 2y + 5. A flat tilted sheet."""
    x, y = point
    return 3.0 * x - 2.0 * y + 5.0


def plane_gradient(point):
    """grad f = (3, -2), at every point in the plane. The point is ignored."""
    del point
    return np.array([3.0, -2.0])


# --------------------------------------------------------------------------
# 3. A product -- the smallest function whose partials involve each other
# --------------------------------------------------------------------------

def product(point):
    """f(x, y) = xy. Flat along both axes through the origin, curved between.

    This is the function that makes the phrase "hold the other one still" do
    real work. Along the x-axis (y = 0) the function is identically zero, so
    the slope in x is zero. Move off that line and the slope in x is y.
    """
    x, y = point
    return x * y


def product_gradient(point):
    """grad f = (y, x). Each partial is the OTHER variable, held fixed."""
    x, y = point
    return np.array([y, x])


# --------------------------------------------------------------------------
# 4. A saddle -- a stationary point that is neither a peak nor a floor
# --------------------------------------------------------------------------

def saddle(point):
    """f(x, y) = x^2 - y^2. Up along x, down along y, flat at the origin."""
    x, y = point
    return x * x - y * y


def saddle_gradient(point):
    """grad f = (2x, -2y). Zero at the origin, and the origin is a saddle."""
    x, y = point
    return np.array([2.0 * x, -2.0 * y])


# --------------------------------------------------------------------------
# 5. A dome -- a stationary point that IS a maximum
# --------------------------------------------------------------------------

def dome(point):
    """f(x, y) = -(x^2 + y^2). The bowl turned upside down."""
    x, y = point
    return -(x * x + y * y)


def dome_gradient(point):
    """grad f = (-2x, -2y). Also zero at the origin. Also a stationary point."""
    x, y = point
    return np.array([-2.0 * x, -2.0 * y])


# --------------------------------------------------------------------------
# 6. A genuine cubic -- the one surface where truncation error is visible
# --------------------------------------------------------------------------

def cubic(point):
    """f(x, y) = x^3 + x*y^2.

    Every other surface here is at most quadratic in each variable, which
    makes the central difference exact for them. This one is not, and its
    error is not merely small but PREDICTABLE: expanding
    ((x+h)^3 - (x-h)^3) / (2h) gives 3x^2 + h^2 exactly, so the numerical
    partial in x overshoots the true one by exactly h squared, with no other
    terms at all. Script 06 measures that and checks it to twelve decimal
    places.
    """
    x, y = point
    return x ** 3 + x * y * y


def cubic_gradient(point):
    """grad f = (3x^2 + y^2, 2xy)."""
    x, y = point
    return np.array([3.0 * x * x + y * y, 2.0 * x * y])


# --------------------------------------------------------------------------
# The registry the scripts and tests iterate over
# --------------------------------------------------------------------------

SURFACES = {
    "bowl": (bowl, bowl_gradient, "x^2 + 3y^2", "grad = (2x, 6y)"),
    "plane": (plane, plane_gradient, "3x - 2y + 5", "grad = (3, -2)"),
    "product": (product, product_gradient, "xy", "grad = (y, x)"),
    "saddle": (saddle, saddle_gradient, "x^2 - y^2", "grad = (2x, -2y)"),
    "dome": (dome, dome_gradient, "-(x^2 + y^2)", "grad = (-2x, -2y)"),
    "cubic": (cubic, cubic_gradient, "x^3 + x*y^2", "grad = (3x^2 + y^2, 2xy)"),
}

# The points every surface is probed at. Chosen by hand: one in the first
# quadrant, one with a negative coordinate, one off-axis with a fraction, one
# ON an axis where a partial vanishes, and one far out.
PROBE_POINTS = (
    (1.0, 1.0),
    (2.0, -1.0),
    (-0.5, 3.0),
    (4.0, 0.0),
    (0.25, 0.75),
)

# The three surfaces that all have a zero gradient at the origin, and what the
# origin actually IS for each. The gradient cannot tell them apart; this table
# is the answer the gradient does not carry.
STATIONARY_AT_ORIGIN = (
    ("bowl", "minimum", "every direction goes up"),
    ("dome", "maximum", "every direction goes down"),
    ("saddle", "saddle", "up along x, down along y"),
)


# --------------------------------------------------------------------------
# Exact contours, parametrised WITHOUT reference to the gradient
# --------------------------------------------------------------------------
#
# The claim being tested is that the gradient is perpendicular to the contour.
# Walking along the contour by stepping perpendicular to the gradient would
# make that claim true by construction and prove nothing. So each contour
# below is an exact algebraic parametrisation, derived on paper from the
# function alone, and the lab checks it lands back on the same value of f
# before it uses it for anything.

def bowl_contour(level, t):
    """A point on the ellipse x^2 + 3y^2 = level, for the angle parameter t.

    Substitute x = sqrt(level) cos t and y = sqrt(level/3) sin t into
    x^2 + 3y^2 and you get level (cos^2 t + sin^2 t) = level, for every t.
    """
    a = np.sqrt(level)
    b = np.sqrt(level / 3.0)
    return np.array([a * np.cos(t), b * np.sin(t)])


def product_contour(level, t):
    """A point on the hyperbola xy = level, for the parameter t (x = t)."""
    return np.array([t, level / t])


def dome_contour(level, t):
    """A point on the circle -(x^2 + y^2) = level, for the angle t.

    level must be negative; the radius is sqrt(-level).
    """
    r = np.sqrt(-level)
    return np.array([r * np.cos(t), r * np.sin(t)])


CONTOURS = {
    "bowl": (bowl, bowl_contour, 4.0, 0.7),
    "product": (product, product_contour, 6.0, 2.0),
    "dome": (dome, dome_contour, -9.0, 1.1),
}


# --------------------------------------------------------------------------
# A three-parameter model, so "one partial per parameter" is not just a claim
# --------------------------------------------------------------------------
#
# Four invented samples. Each row is (a, b, target). Nothing was measured to
# produce these; they were chosen so the arithmetic stays checkable by hand.

SAMPLES = (
    (1.0, 2.0, 8.0),
    (2.0, 1.0, 7.0),
    (3.0, 3.0, 15.0),
    (0.0, 1.0, 3.0),
)

# The parameter vector the lab evaluates the loss and its gradient at.
START_PARAMS = (1.0, 1.0, 1.0)


def model_loss(params):
    """Mean squared error of `pred = w1*a + w2*b + c` over SAMPLES.

    Three parameters, so the gradient is a vector of three numbers -- one per
    parameter. A real network has millions of parameters and the gradient has
    millions of entries. The idea does not change; only the count does.
    """
    w1, w2, c = params
    total = 0.0
    for a, b, target in SAMPLES:
        residual = w1 * a + w2 * b + c - target
        total += residual * residual
    return total / len(SAMPLES)


def model_loss_gradient(params):
    """The exact gradient, differentiated by hand.

    L = (1/n) sum r_i^2 with r_i = w1*a_i + w2*b_i + c - y_i, so
    dL/dw1 = (2/n) sum r_i * a_i, dL/dw2 = (2/n) sum r_i * b_i,
    dL/dc  = (2/n) sum r_i.
    """
    w1, w2, c = params
    n = len(SAMPLES)
    g_w1 = g_w2 = g_c = 0.0
    for a, b, target in SAMPLES:
        residual = w1 * a + w2 * b + c - target
        g_w1 += 2.0 * residual * a
        g_w2 += 2.0 * residual * b
        g_c += 2.0 * residual
    return np.array([g_w1 / n, g_w2 / n, g_c / n])

"""The reference suite: every claim this lab makes, checked against a value.

Run from the lab directory:

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

Nothing here reads source code or checks that a file exists. Every test calls
something and compares the answer to a number that was worked out by hand or
derived algebraically, with a tolerance stated in `surfaces.py` and justified
there rather than tuned until green.
"""

from __future__ import annotations

import math

import numpy as np
import pytest

import surfaces as S
from gradients import (
    angle_degrees,
    angular_gap_degrees,
    contour_chord,
    directional_derivative,
    directional_derivative_direct,
    forward_partial,
    gradient,
    magnitude,
    partial,
    sweep_directions,
    unit,
)

ALL_SURFACES = sorted(S.SURFACES)


# ==========================================================================
# 1. partial: one input moves, the rest hold still
# ==========================================================================

@pytest.mark.parametrize("name", ALL_SURFACES)
@pytest.mark.parametrize("point", S.PROBE_POINTS)
@pytest.mark.parametrize("index", (0, 1))
def test_partial_matches_the_hand_derived_value(name, point, index):
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    assert partial(f, point, index) == pytest.approx(
        exact_gradient(point)[index], abs=S.GRADIENT_TOL
    )


def test_partial_leaves_the_other_coordinates_untouched():
    """The defining property. If the y coordinate moved, this is not a partial."""
    seen = []

    def spy(p):
        seen.append(tuple(float(v) for v in p))
        return S.product(p)

    partial(spy, (2.0, 5.0), 0)
    assert len(seen) == 2
    assert {p[1] for p in seen} == {5.0}, "y must be identical in both evaluations"
    assert sorted(p[0] for p in seen) == [2.0 - S.H_DEFAULT, 2.0 + S.H_DEFAULT]


def test_partial_does_not_mutate_the_point_it_was_given():
    point = np.array([1.0, 2.0])
    partial(S.bowl, point, 0)
    assert point.tolist() == [1.0, 2.0]


def test_partial_accepts_a_tuple_a_list_and_an_array_alike():
    answers = [
        partial(S.bowl, (1.0, 1.0), 0),
        partial(S.bowl, [1.0, 1.0], 0),
        partial(S.bowl, np.array([1.0, 1.0]), 0),
    ]
    assert answers[0] == answers[1] == answers[2]


def test_central_difference_is_exact_on_a_quadratic_for_any_step():
    """((x+h)^2 - (x-h)^2) / 2h = 2x algebraically, so h barely matters."""
    for h in (1e-1, 1e-2, 1e-3, 1e-4, 1e-5):
        assert partial(S.bowl, (2.0, 1.0), 0, h) == pytest.approx(4.0, abs=1e-7)


def test_a_zero_partial_does_not_mean_a_flat_point():
    """At (1, 0) on f = xy the x-slope is zero and the y-slope is 1."""
    assert partial(S.product, (1.0, 0.0), 0) == pytest.approx(0.0, abs=S.GRADIENT_TOL)
    assert partial(S.product, (1.0, 0.0), 1) == pytest.approx(1.0, abs=S.GRADIENT_TOL)


# ==========================================================================
# 2. gradient: the vector of partials
# ==========================================================================

@pytest.mark.parametrize("name", ALL_SURFACES)
@pytest.mark.parametrize("point", S.PROBE_POINTS)
def test_gradient_matches_the_hand_derived_vector(name, point):
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    assert gradient(f, point) == pytest.approx(
        exact_gradient(point), abs=S.GRADIENT_TOL
    )


@pytest.mark.parametrize("name", ALL_SURFACES)
@pytest.mark.parametrize("point", S.PROBE_POINTS)
def test_gradient_is_exactly_the_partials_side_by_side(name, point):
    f = S.SURFACES[name][0]
    built = np.array([partial(f, point, 0), partial(f, point, 1)])
    assert gradient(f, point) == pytest.approx(built, abs=0.0)


def test_gradient_has_one_entry_per_input_not_one_per_dimension_of_the_surface():
    assert gradient(S.bowl, (1.0, 1.0)).size == 2
    assert gradient(S.model_loss, (1.0, 1.0, 1.0)).size == 3


def test_gradient_works_on_three_inputs():
    assert gradient(S.model_loss, S.START_PARAMS) == pytest.approx(
        S.model_loss_gradient(S.START_PARAMS), abs=S.GRADIENT_TOL
    )


def test_the_model_gradient_is_the_three_whole_numbers_worked_out_by_hand():
    assert S.model_loss(S.START_PARAMS) == 22.5
    assert S.model_loss_gradient(S.START_PARAMS).tolist() == [-17.0, -18.0, -8.0]


def test_gradient_returns_a_numpy_array_of_floats():
    g = gradient(S.bowl, (1.0, 1.0))
    assert isinstance(g, np.ndarray)
    assert g.dtype == np.float64


# ==========================================================================
# 3. magnitude and unit
# ==========================================================================

def test_magnitude_is_the_euclidean_norm_from_day_99():
    assert magnitude([3.0, 4.0]) == 5.0
    assert magnitude([1.0, 1.0]) == pytest.approx(math.sqrt(2.0))
    assert magnitude([0.0, 0.0]) == 0.0


def test_magnitude_of_the_bowl_gradient_at_one_one():
    """grad = (2, 6), so the length is sqrt(4 + 36) = sqrt(40)."""
    assert magnitude(gradient(S.bowl, (1.0, 1.0))) == pytest.approx(
        math.sqrt(40.0), abs=1e-9
    )


@pytest.mark.parametrize("point", S.PROBE_POINTS[:3])
def test_unit_has_length_one_and_keeps_the_bearing(point):
    g = gradient(S.cubic, point)
    u = unit(g)
    assert magnitude(u) == pytest.approx(1.0, abs=1e-12)
    assert angle_degrees(u) == pytest.approx(angle_degrees(g), abs=1e-9)


def test_unit_refuses_the_zero_vector_rather_than_dividing_by_zero():
    with pytest.raises(ValueError, match="no direction"):
        unit([0.0, 0.0])


def test_unit_of_a_long_arrow_and_a_short_one_in_the_same_direction_agree():
    assert unit([1.0, 2.0]) == pytest.approx(unit([50.0, 100.0]), abs=1e-15)


# ==========================================================================
# 4. directional derivatives -- Day 103's dot product doing real work
# ==========================================================================

DIRECTIONS = ((1.0, 0.0), (0.0, 1.0), (1.0, 1.0), (-1.0, 2.0),
              (3.0, -1.0), (-2.0, -5.0), (7.0, 0.5))


@pytest.mark.parametrize("direction", DIRECTIONS)
@pytest.mark.parametrize("name", ("bowl", "product", "cubic", "saddle"))
def test_dot_with_the_gradient_agrees_with_a_direct_measurement(name, direction):
    """The claim that grad . u IS the rate of change, checked without assuming it."""
    f = S.SURFACES[name][0]
    point = (1.0, 1.0)
    assert directional_derivative(f, point, direction) == pytest.approx(
        directional_derivative_direct(f, point, direction), abs=S.GRADIENT_TOL
    )


def test_a_partial_derivative_is_the_directional_derivative_along_an_axis():
    point = (1.0, 1.0)
    assert directional_derivative(S.bowl, point, (1.0, 0.0)) == pytest.approx(
        2.0, abs=S.GRADIENT_TOL
    )
    assert directional_derivative(S.bowl, point, (0.0, 1.0)) == pytest.approx(
        6.0, abs=S.GRADIENT_TOL
    )


def test_the_length_of_the_direction_arrow_does_not_change_the_answer():
    point = (1.0, 1.0)
    short = directional_derivative(S.bowl, point, (1.0, 2.0))
    long = directional_derivative(S.bowl, point, (1000.0, 2000.0))
    assert short == pytest.approx(long, abs=1e-9)


def test_walking_backwards_negates_the_rate():
    point = (2.0, -1.0)
    forward = directional_derivative(S.cubic, point, (1.0, 3.0))
    backward = directional_derivative(S.cubic, point, (-1.0, -3.0))
    assert forward == pytest.approx(-backward, abs=1e-9)


def test_a_direction_perpendicular_to_the_gradient_gives_zero():
    """grad of the bowl at (1, 1) is (2, 6); (3, -1) dots to 3*2 - 1*6 = 0."""
    assert directional_derivative(S.bowl, (1.0, 1.0), (3.0, -1.0)) == pytest.approx(
        0.0, abs=S.GRADIENT_TOL
    )


# ==========================================================================
# 5. steepest ascent -- the first fact that must be demonstrated
# ==========================================================================

SWEEP_TRIALS = (
    ("bowl", (1.0, 1.0)),
    ("bowl", (0.25, 0.75)),
    ("bowl", (3.0, -2.0)),
    ("bowl", (1.0, 0.4)),
    ("product", (2.0, -1.0)),
    ("saddle", (1.5, 0.5)),
    ("cubic", (1.0, 1.0)),
    ("plane", (-2.0, 4.0)),
)


@pytest.mark.parametrize("name,point", SWEEP_TRIALS)
def test_the_best_of_360_directions_is_the_gradient_direction(name, point):
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    angles, rates = sweep_directions(f, point)
    best = int(np.argmax(rates))
    gap = angular_gap_degrees(
        float(np.degrees(angles[best])), angle_degrees(exact_gradient(point))
    )
    assert gap <= S.ANGLE_TOL_DEGREES


@pytest.mark.parametrize("name,point", SWEEP_TRIALS)
def test_no_direction_at_all_beats_the_gradients_own_magnitude(name, point):
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    _angles, rates = sweep_directions(f, point)
    steepest = magnitude(exact_gradient(point))
    assert float(np.max(rates)) <= steepest + S.GRADIENT_TOL


@pytest.mark.parametrize("name,point", SWEEP_TRIALS)
def test_the_winning_rate_is_the_magnitude_times_cosine_of_the_sampling_gap(name, point):
    """Day 103's geometric dot product, to nine decimal places."""
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    angles, rates = sweep_directions(f, point)
    best = int(np.argmax(rates))
    gap = angular_gap_degrees(
        float(np.degrees(angles[best])), angle_degrees(exact_gradient(point))
    )
    ratio = float(np.max(rates)) / magnitude(exact_gradient(point))
    assert ratio == pytest.approx(math.cos(math.radians(gap)), abs=1e-9)


def test_the_steepest_descent_is_the_exact_opposite_bearing():
    angles, rates = sweep_directions(S.bowl, (1.0, 1.0))
    up = float(np.degrees(angles[int(np.argmax(rates))]))
    down = float(np.degrees(angles[int(np.argmin(rates))]))
    assert angular_gap_degrees(up, down) == pytest.approx(180.0, abs=1e-9)
    assert float(np.max(rates)) == pytest.approx(-float(np.min(rates)), abs=1e-6)


def test_the_two_flattest_bearings_are_ninety_degrees_from_the_gradient():
    point = (1.0, 1.0)
    angles, rates = sweep_directions(S.bowl, point)
    bearing = angle_degrees(S.bowl_gradient(point))
    for i in np.argsort(np.abs(rates))[:2]:
        off = angular_gap_degrees(float(np.degrees(angles[i])), bearing)
        assert off == pytest.approx(90.0, abs=S.ANGLE_TOL_DEGREES)


def test_the_residual_gap_is_bounded_by_the_sampling_not_by_the_calculus():
    """The gap can never exceed half the spacing between sampled bearings.

    Note what this test does NOT assert. The obvious version -- that a finer
    sweep always gives a strictly smaller gap -- is false, and was found to be
    false here: the bowl's gradient at (1, 1) has bearing 71.5651 degrees, and
    both a 60-direction sweep (every 6 degrees) and a 360-direction sweep
    (every 1 degree) land on 72, so both leave exactly the same 0.4349-degree
    gap. Sampling more finely guarantees a smaller BOUND, not a smaller gap on
    any particular bearing.
    """
    point = (1.0, 1.0)
    bearing = angle_degrees(S.bowl_gradient(point))
    gaps = {}
    for n in (60, 360, 3600):
        angles, rates = sweep_directions(S.bowl, point, n=n)
        best = int(np.argmax(rates))
        gap = angular_gap_degrees(float(np.degrees(angles[best])), bearing)
        gaps[n] = gap
        assert gap <= 180.0 / n, (n, gap)
    assert gaps[60] == pytest.approx(gaps[360], abs=1e-12)
    assert gaps[3600] < gaps[360]
    assert gaps[3600] < 0.05


# ==========================================================================
# 6. perpendicular to the contour -- the second fact
# ==========================================================================

@pytest.mark.parametrize("name", sorted(S.CONTOURS))
@pytest.mark.parametrize("t", (0.4, 0.9, 1.4, 1.9))
def test_the_parametrised_contour_really_does_hold_f_constant(name, t):
    """Checked BEFORE the contour is used, so the algebra is not taken on trust."""
    f, contour, level, _t0 = S.CONTOURS[name]
    assert f(contour(level, t)) == pytest.approx(level, abs=1e-12)


@pytest.mark.parametrize("name", sorted(S.CONTOURS))
@pytest.mark.parametrize("t", (0.4, 0.9, 1.4, 1.9))
def test_the_gradient_is_perpendicular_to_a_step_along_the_contour(name, t):
    f, contour, level, _t0 = S.CONTOURS[name]
    chord, p, _q, _fp, _fq = contour_chord(f, contour, level, t, S.CONTOUR_DELTA)
    assert abs(float(np.dot(unit(gradient(f, p)), chord))) < S.CONTOUR_DOT_TOL


@pytest.mark.parametrize("name", sorted(S.CONTOURS))
def test_the_dot_product_shrinks_in_step_with_the_contour_step(name):
    """The honest form of 'it goes to zero': first-order, ten for ten."""
    f, contour, level, t0 = S.CONTOURS[name]
    previous = None
    for k in (2, 3, 4, 5, 6):
        delta = 10.0 ** (-k)
        chord, p, _q, _fp, _fq = contour_chord(f, contour, level, t0, delta)
        dot = abs(float(np.dot(unit(gradient(f, p)), chord)))
        if previous is not None:
            assert 9.0 < previous / dot < 11.0
        previous = dot


def test_the_exact_tangent_and_the_exact_gradient_dot_to_exactly_zero():
    """No tolerance needed: the algebra cancels term for term."""
    level = 4.0
    a = math.sqrt(level)
    b = math.sqrt(level / 3.0)
    for t in (0.0, 0.4, 0.9, 1.4, 1.9, 2.7):
        p = S.bowl_contour(level, t)
        tangent = np.array([-a * math.sin(t), b * math.cos(t)])
        assert abs(float(np.dot(tangent, S.bowl_gradient(p)))) < 1e-14


def test_a_step_along_the_contour_changes_f_far_less_than_a_step_across_it():
    point = np.array([1.0, 1.0])
    g = gradient(S.bowl, point)
    along = unit(np.array([-g[1], g[0]]))
    across = unit(g)
    step = 0.001
    gain_along = abs(S.bowl(point + step * along) - S.bowl(point))
    gain_across = abs(S.bowl(point + step * across) - S.bowl(point))
    assert gain_across > 100 * gain_along


def test_the_gradients_length_predicts_the_gain_of_a_small_step_up_it():
    point = np.array([1.0, 1.0])
    g = gradient(S.bowl, point)
    step = 0.001
    measured = S.bowl(point + step * unit(g)) - S.bowl(point)
    assert measured == pytest.approx(magnitude(g) * step, abs=1e-4)


# ==========================================================================
# 7. linear and quadratic
# ==========================================================================

@pytest.mark.parametrize("point", ((0.0, 0.0), (1.0, 1.0), (-40.0, 17.5),
                                   (0.001, 0.002), (7.5, -3.25)))
def test_the_gradient_of_a_plane_is_the_same_vector_everywhere(point):
    assert gradient(S.plane, point) == pytest.approx([3.0, -2.0], abs=S.GRADIENT_TOL)


def test_the_planes_gradient_never_varies_between_points():
    seen = [gradient(S.plane, p) for p in ((0.0, 0.0), (1.0, 1.0), (-40.0, 17.5))]
    assert float(np.max(np.abs(np.array(seen) - seen[0]))) < S.GRADIENT_TOL


@pytest.mark.parametrize("point", ((0.5, 0.5), (1.0, 1.0), (2.0, 2.0),
                                   (4.0, 4.0), (-3.0, 1.0)))
def test_the_bowls_gradient_points_away_from_its_minimum(point):
    outward = float(np.dot(unit(gradient(S.bowl, point)), unit(np.array(point))))
    assert outward > 0.0


def test_the_bowls_gradient_grows_with_distance_from_the_minimum():
    lengths = [magnitude(gradient(S.bowl, (r, r))) for r in (0.5, 1.0, 2.0, 4.0)]
    assert lengths == sorted(lengths)
    assert lengths[1] == pytest.approx(2.0 * lengths[0], rel=1e-6)


def test_the_negative_gradient_does_not_point_straight_at_the_minimum_on_an_ellipse():
    """The mismatch that makes gradient descent zig-zag. Real, and measured."""
    point = (1.0, 1.0)
    back = angle_degrees(-gradient(S.bowl, point))
    straight = angle_degrees(-np.array(point))
    assert angular_gap_degrees(back, straight) > 20.0


# ==========================================================================
# 8. the zero gradient, and what it cannot tell you
# ==========================================================================

@pytest.mark.parametrize("name,kind,_why", S.STATIONARY_AT_ORIGIN)
def test_all_three_surfaces_have_a_zero_gradient_at_the_origin(name, kind, _why):
    del kind, _why
    assert magnitude(gradient(S.SURFACES[name][0], (0.0, 0.0))) < S.GRADIENT_TOL


def test_the_three_zero_gradients_are_indistinguishable_from_each_other():
    vectors = [gradient(S.SURFACES[n][0], (0.0, 0.0))
               for n, _k, _w in S.STATIONARY_AT_ORIGIN]
    for v in vectors[1:]:
        assert v == pytest.approx(vectors[0], abs=S.GRADIENT_TOL)


def test_but_the_points_themselves_are_completely_different():
    radius = 0.1
    verdicts = {}
    for name, _kind, _why in S.STATIONARY_AT_ORIGIN:
        f = S.SURFACES[name][0]
        changes = [f(radius * np.array([math.cos(a), math.sin(a)])) - f(np.zeros(2))
                   for a in np.radians(np.arange(0, 360, 45))]
        up = sum(1 for c in changes if c > 1e-12)
        down = sum(1 for c in changes if c < -1e-12)
        verdicts[name] = (up > 0, down > 0)
    assert verdicts["bowl"] == (True, False)     # a minimum: only up
    assert verdicts["dome"] == (False, True)     # a maximum: only down
    assert verdicts["saddle"] == (True, True)    # a saddle: both


def test_the_saddle_rises_along_x_and_falls_along_y():
    assert S.saddle((0.5, 0.0)) > 0.0
    assert S.saddle((0.0, 0.5)) < 0.0
    assert S.saddle((0.5, 0.5)) == 0.0


def test_near_a_saddle_the_gradient_is_small_without_being_zero():
    lengths = [magnitude(gradient(S.saddle, (r, r))) for r in (1.0, 0.1, 0.01)]
    assert lengths[0] > lengths[1] > lengths[2] > 0.0


# ==========================================================================
# 9. step size, and Day 108's U-curve
# ==========================================================================

@pytest.mark.parametrize("h", (1e-1, 1e-2, 1e-3))
def test_the_cubics_truncation_error_is_exactly_h_squared(h):
    exact = float(S.cubic_gradient((2.0, 1.0))[0])
    error = partial(S.cubic, (2.0, 1.0), 0, h) - exact
    assert error == pytest.approx(h * h, rel=1e-5)


def test_the_error_curve_is_u_shaped_rather_than_monotonic():
    exact = float(S.cubic_gradient((2.0, 1.0))[0])
    errors = {k: abs(partial(S.cubic, (2.0, 1.0), 0, 10.0 ** -k) - exact)
              for k in range(0, 15)}
    best = min(errors, key=errors.get)
    assert 0 < best < 14, "the best step is in the middle, not at either end"
    assert errors[0] > errors[best] < errors[14]


def test_the_default_step_sits_at_the_bottom_of_that_curve():
    exact = float(S.cubic_gradient((2.0, 1.0))[0])
    errors = {k: abs(partial(S.cubic, (2.0, 1.0), 0, 10.0 ** -k) - exact)
              for k in range(0, 15)}
    assert 10.0 ** -min(errors, key=errors.get) == S.H_DEFAULT


def test_too_small_a_step_is_worse_than_a_step_a_trillion_times_bigger():
    exact = float(S.cubic_gradient((2.0, 1.0))[0])
    tiny = abs(partial(S.cubic, (2.0, 1.0), 0, 1e-14) - exact)
    large = abs(partial(S.cubic, (2.0, 1.0), 0, 1e-1) - exact)
    assert tiny > large


def test_the_central_difference_beats_the_forward_one_at_the_default_step():
    exact = float(S.cubic_gradient((2.0, 1.0))[0])
    central = abs(partial(S.cubic, (2.0, 1.0), 0, S.H_DEFAULT) - exact)
    forward = abs(forward_partial(S.cubic, (2.0, 1.0), 0, S.H_DEFAULT) - exact)
    assert forward > 1000 * central


def test_the_forward_differences_error_falls_like_h_not_h_squared():
    exact = float(S.cubic_gradient((2.0, 1.0))[0])
    errors = [abs(forward_partial(S.cubic, (2.0, 1.0), 0, h) - exact)
              for h in (1e-1, 1e-2, 1e-3)]
    assert errors[0] / errors[1] == pytest.approx(10.0, rel=0.05)
    assert errors[1] / errors[2] == pytest.approx(10.0, rel=0.05)


def test_roundoff_grows_with_the_size_of_f_and_the_bound_predicts_it():
    """Found by watching an assertion fail at (1000, -1000). Kept, and checked."""
    eps = float(np.finfo(float).eps)
    for point in ((1000.0, -1000.0), (100000.0, -100000.0), (10000000.0, -10000000.0)):
        error = abs(partial(S.plane, point, 0) - 3.0)
        predicted = eps * abs(S.plane(point)) / (2.0 * S.H_DEFAULT)
        assert error < 3.0 * predicted
        assert error > predicted / 100.0


def test_the_numerical_gradient_stops_meeting_the_labs_tolerance_far_from_home():
    """The honest boundary on every other tolerance in this file."""
    assert abs(partial(S.plane, (1000.0, -1000.0), 0) - 3.0) > S.GRADIENT_TOL


# ==========================================================================
# 10. numpy.gradient does a related, different job
# ==========================================================================

def _grid():
    xs = np.linspace(0.0, 4.0, 9)
    ys = np.linspace(0.0, 4.0, 9)
    gx, gy = np.meshgrid(xs, ys, indexing="ij")
    return xs, ys, gx, gy


def test_numpy_gradient_is_exact_in_the_interior_of_a_sampled_quadratic():
    xs, ys, X, Y = _grid()
    gx, gy = np.gradient(X * X + 3.0 * Y * Y, xs, ys)
    assert gx[2, 2] == pytest.approx(2.0 * xs[2], abs=1e-12)
    assert gy[2, 2] == pytest.approx(6.0 * ys[2], abs=1e-12)


def test_numpy_gradient_defaults_to_first_order_at_the_boundary():
    """A real default worth knowing: the corner of an exact quadratic is wrong."""
    xs, ys, X, Y = _grid()
    gx, gy = np.gradient(X * X + 3.0 * Y * Y, xs, ys)
    assert gx[0, 0] == pytest.approx(0.5, abs=1e-12)
    assert gy[0, 0] == pytest.approx(1.5, abs=1e-12)
    assert gx[0, 0] != pytest.approx(0.0, abs=1e-6)


def test_edge_order_two_fixes_that_corner_exactly():
    xs, ys, X, Y = _grid()
    gx, gy = np.gradient(X * X + 3.0 * Y * Y, xs, ys, edge_order=2)
    assert gx[0, 0] == pytest.approx(0.0, abs=1e-12)
    assert gy[0, 0] == pytest.approx(0.0, abs=1e-12)


def test_numpy_gradients_error_on_a_cubic_is_the_grid_spacing_squared():
    """Same h^2 law -- but h is the sampling, which you cannot choose."""
    xs, ys, X, Y = _grid()
    spacing = float(xs[1] - xs[0])
    gx, _gy = np.gradient(X ** 3 + X * Y * Y, xs, ys, edge_order=2)
    exact = 3.0 * xs[4] ** 2 + ys[4] ** 2
    assert abs(gx[4, 4] - exact) == pytest.approx(spacing ** 2, abs=1e-12)


def test_our_gradient_gets_the_same_point_right_because_it_chooses_its_own_step():
    xs, ys, _X, _Y = _grid()
    point = (float(xs[4]), float(ys[4]))
    assert gradient(S.cubic, point) == pytest.approx(
        S.cubic_gradient(point), abs=S.GRADIENT_TOL
    )


def test_numpy_gradient_returns_a_field_not_a_vector():
    """The shape difference IS the conceptual difference."""
    xs, ys, X, Y = _grid()
    gx, _gy = np.gradient(X * X + 3.0 * Y * Y, xs, ys)
    assert gx.shape == (9, 9)
    assert gradient(S.bowl, (1.0, 1.0)).shape == (2,)


# ==========================================================================
# 11. one partial per parameter
# ==========================================================================

def test_the_loss_is_the_mean_of_the_four_squared_residuals():
    w1, w2, c = S.START_PARAMS
    by_hand = sum((w1 * a + w2 * b + c - y) ** 2 for a, b, y in S.SAMPLES) / 4
    assert S.model_loss(S.START_PARAMS) == by_hand == 22.5


def test_a_numerical_gradient_costs_exactly_two_evaluations_per_parameter():
    calls = []

    def counted(p):
        calls.append(1)
        return S.model_loss(p)

    gradient(counted, S.START_PARAMS)
    assert len(calls) == 2 * len(S.START_PARAMS)


def test_a_small_step_against_the_gradient_reduces_the_loss():
    before = S.model_loss(S.START_PARAMS)
    moved = np.array(S.START_PARAMS) - 0.01 * S.model_loss_gradient(S.START_PARAMS)
    assert S.model_loss(moved) < before


def test_a_small_step_ALONG_the_gradient_increases_it():
    before = S.model_loss(S.START_PARAMS)
    moved = np.array(S.START_PARAMS) + 0.01 * S.model_loss_gradient(S.START_PARAMS)
    assert S.model_loss(moved) > before


def test_too_large_a_step_overshoots_and_ends_up_worse_than_the_start():
    before = S.model_loss(S.START_PARAMS)
    moved = np.array(S.START_PARAMS) - 0.2 * S.model_loss_gradient(S.START_PARAMS)
    assert S.model_loss(moved) > before


def test_the_three_parameter_gradient_still_obeys_steepest_ascent():
    """No picture, same fact: no direction beats the gradient's own length."""
    rng = np.random.default_rng(S.SEED)
    point = S.START_PARAMS
    g = S.model_loss_gradient(point)
    best = magnitude(g)
    for _ in range(200):
        direction = unit(rng.normal(size=3))
        rate = directional_derivative_direct(S.model_loss, point, direction)
        assert rate <= best + 1e-6


def test_the_seeded_generator_is_stated_and_reproducible():
    first = np.random.default_rng(S.SEED).normal(size=3)
    second = np.random.default_rng(S.SEED).normal(size=3)
    assert first.tolist() == second.tolist()


# ==========================================================================
# 12. housekeeping
# ==========================================================================

def test_every_surface_in_the_registry_has_a_matching_exact_gradient():
    for name, (f, exact_gradient, expression, gradient_expression) in S.SURFACES.items():
        assert callable(f) and callable(exact_gradient), name
        assert expression and gradient_expression, name
        assert gradient(f, (1.0, 1.0)) == pytest.approx(
            exact_gradient((1.0, 1.0)), abs=S.GRADIENT_TOL
        ), name


def test_the_tolerances_are_stated_in_one_place_and_are_not_zero():
    for value in (S.H_DEFAULT, S.GRADIENT_TOL, S.CONTOUR_DELTA,
                  S.CONTOUR_DOT_TOL, S.ANGLE_TOL_DEGREES):
        assert value > 0.0


def test_the_gradient_tolerance_has_real_headroom_over_the_worst_probe():
    worst = 0.0
    for name in ALL_SURFACES:
        f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
        for point in S.PROBE_POINTS:
            worst = max(worst, float(np.max(np.abs(
                gradient(f, point) - exact_gradient(point)))))
    assert worst < S.GRADIENT_TOL
    assert S.GRADIENT_TOL / worst > 10.0, "the tolerance is not merely scraping past"

lesson_id: D109
day: 109
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-109-partial-derivatives-and-gradients
  - 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_hold_everything_else_still.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_the_gradient_vector.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_steepest_ascent.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_perpendicular_to_the_contour.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_flat_ground_three_ways.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_step_size_and_the_u_curve.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 07_one_partial_per_parameter.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 . -path ./.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 -> 271 passed; pytest starter -> 1 passed, 205 skipped on an untouched checkout, and 206 passed against a fully solved copy of starter/ kept outside the lab. All seven reference scripts exit 0 with every internal assertion holding. Everything was run through a real lab-local .venv created by the documented setup commands, not through an authoring environment. Network is needed once to install numpy and pytest; nothing else in the lab opens a socket, and section 7 of the harness greps the sources to prove it. Section 6 re-runs the harness with one expectation deliberately swapped for the belief that a bowl''s gradient at (1, 1) points at 45 degrees -- straight away from the minimum -- rather than at its true 71.5651, 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 .venv defect fixed across six labs earlier today was verified NOT reintroduced: the whole suite was re-run in a copy of the lab with no .venv at all (exit 0, 98 checks) and again in that copy with a fake .venv containing foreign __pycache__ and .pytest_cache directories (exit 0, 98 checks), and an unpruned find was confirmed to see those directories while the pruned find in section 7 does not; .venv is never reported as a stray file. Measured on this run: thirty numerical gradients across six surfaces and five points all agree with hand-derived exact gradients, worst single error 3.961e-10 against an asserted tolerance of 1e-8 that the suite also requires to have at least tenfold headroom. Of 360 bearings measured by direct central difference with no gradient involved, the fastest climb is at 72 degrees against the gradient''s true bearing of 71.5651, a gap of 0.4349 which is inside the half-degree the sampling allows, and the winning rate divided by the gradient magnitude equals cos(gap) to better than 1e-9. Perpendicularity to three exactly parametrised contours gives a worst dot product of 6.645e-06 against a 1e-4 tolerance, and the dot product divides by ten each time the contour step divides by ten (measured ratios 9.9758 to 10.0013). Three unplanned findings are kept rather than tidied away: a plane''s gradient, which is constant, cannot be measured to the lab''s own tolerance at (1000, -1000), and the roundoff bound eps*|f|/(2h) tracks the measured error across seven orders of magnitude, so the suite asserts that failure instead of widening the tolerance; a 60-direction and a 360-direction sweep leave EXACTLY the same 0.4349-degree gap at that point because both grids contain 72, so the test asserts the real law (the gap never exceeds half the sampling step) rather than the plausible but false claim that a finer sweep always does better; and numpy.gradient with its default edge_order=1 returns (0.5, 1.5) at the corner of an exactly sampled quadratic where (0, 0) is correct, which is asserted so a future release changing that default would fail the suite rather than let this page go stale.'

# What this lab installs, and what it costs you

Two packages, both free, both open source, no account and no key.

| Package | Version pinned here | Licence | Why it is here |
| --- | --- | --- | --- |
| `numpy` | 2.5.2 | BSD 3-Clause | Vectors, dot products, `linspace`, and the trigonometry used to sweep 360 directions around a circle. The gradient machinery itself is written from scratch; NumPy holds the results. |
| `pytest` | 9.1.1 | MIT | Runs both suites: the reference tests in `examples/` and your running score in `starter/`. |

Both versions are pinned exactly. Section 1 of `tests/run_tests.sh` reads the
installed numpy and compares it against this file rather than trusting it, so a
mismatch is reported rather than discovered later as a puzzling number.

## What NumPy is and is not doing here

It is worth being precise, because NumPy has a function called `gradient` and
this lab has one too.

**NumPy is not computing any derivative in this lab.** Every partial
derivative, every gradient and every directional derivative comes from
`examples/gradients.py`, which evaluates the function at two points and
subtracts. NumPy supplies arrays, `np.dot`, `np.linspace`, `np.cos` and
`np.sqrt`, and holds the answers.

`numpy.gradient` appears once, in script 06, and it is there to be
*distinguished* from ours rather than used. It differences an array of values
already sampled on a grid; ours differences a function it can call at any point
it likes. Section 4 of that script measures the consequence: on a cubic,
NumPy's error is the grid spacing squared, which is fixed by the data you were
given, while ours picks its own step and lands ten orders of magnitude closer.

Two facts about `numpy.gradient` are asserted rather than described, so that a
future release changing either would fail the suite instead of quietly making
this page wrong:

- Interior values on a sampled quadratic come out **exact**, because a central
  difference is algebraically exact for a quadratic at any spacing.
- Boundary values default to a **first-order** one-sided formula, so the corner
  of that same exact quadratic comes out as `(0.5, 1.5)` where `(0, 0)` is
  correct. Passing `edge_order=2` fixes it exactly. That default costs people
  an afternoon reasonably often.

## The network

Installing these two packages is the only thing in this lab that touches the
network. Nothing here opens a socket, reads a URL or needs an API key, and
section 7 of the test harness greps every source file in `examples/` and
`starter/` to prove it.

## If you cannot install anything at all

You can do more of this lab than you might expect, and it would be dishonest to
pretend you can do all of it.

What works on a bare `python3` with only the standard library, if you replace
the handful of `np` calls with `math` and plain lists:

- `partial` and `forward_partial`, which are two evaluations and a subtraction;
- `gradient`, which is a loop over `partial`;
- `magnitude` and `unit`, which need `math.sqrt`;
- every prediction in `starter/answers.py`, which is where most of the thinking
  lives and none of which requires running anything.

What you lose is the sweep of 360 directions, the contour work, and the
comparison with `numpy.gradient` — that is, most of the *evidence*, though not
most of the *reasoning*.

## Disk

Roughly 60 MB for the virtual environment, almost all of it NumPy. `rm -rf
.venv` from the lab directory is a complete undo.

numpy==2.5.2
pytest==9.1.1

# Which Way Is Uphill? — the eight exercises, in order

Work through these in order. Check yourself at any point from the **lab
directory** (the one above this file):

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

On an untouched checkout that prints `1 passed, 205 skipped`. A skip means
"not attempted"; a failure means "attempted and wrong", and prints both your
answer and the right one. When it prints `206 passed`, you are finished.

Do not open `examples/` until you have tried. It is the answer key, and it is
worth reading afterwards — every script there prints its working — but reading
it first turns a lab into a tutorial.

---

## Exercise 1 — the machinery (`starter/gradients.py`)

Eight functions. The first two are the day; the rest are built out of them.

| # | Function | What it does |
| --- | --- | --- |
| 1.1 | `partial(f, point, index, h)` | Nudge ONE coordinate up and down by `h`, hold the rest still, divide by `2h`. |
| 1.2 | `gradient(f, point, h)` | One partial per input, collected into a vector. |
| 1.3 | `magnitude(vector)` | Day 99's Euclidean norm. |
| 1.4 | `unit(vector)` | The same direction at length 1; `ValueError` on the zero vector. |
| 1.5 | `directional_derivative(f, point, direction, h)` | Dot the gradient with the unit direction. |
| 1.6 | `directional_derivative_direct(...)` | The same number, measured along the direction, with no gradient involved. |
| 1.7 | `sweep_directions(f, point, n, h)` | `n` bearings round the circle, each one's rate of change. |
| 1.8 | `forward_partial(...)` | The one-sided version, for the comparison in exercise 6. |

Write 1.1 first and get its tests green before anything else. Everything below
it inherits its errors, and three of its tests exist to catch the three
specific ways it goes wrong: dividing by `h` instead of `2h`, letting a second
coordinate move, and mutating the caller's point.

1.5 and 1.6 compute the same quantity by completely different routes. The test
that compares them is the most important one in the lab: it is the evidence
that dotting the gradient with a direction really does answer "how fast does
`f` change if I walk that way", rather than being a rule you were told.

Once all eight are written, a further thirty-two tests run automatically. They are
the two facts the day exists for — that the gradient wins a sweep of 360
directions, and that it is perpendicular to a contour — plus the plane's
constant gradient, the bowl's outward-pointing one, the three zero gradients,
the `h^2` law and the U-curve. You do not write anything more for those; they
use the functions you already wrote.

## Exercises 2 to 8 — the predictions (`starter/answers.py`)

Fifty-one predictions. Replace each `None`. Almost every one can be done on
paper, and the ones that cannot are asking you to reason about a shape rather
than compute a number.

| # | Topic | What you are asked for |
| --- | --- | --- |
| 2 | Partial derivatives by hand | Eight: `df/dx` and `df/dy` for three surfaces, plus what the rounded `d` is for. |
| 3 | The gradient as a vector | Six: assembling it, its length, how many components it has, and what the length means. |
| 4 | Directional derivatives | Seven: rates along given bearings, the largest and smallest available, and which trigonometric function relates them. |
| 5 | Contours | Six: the shape of the contours, the angle to the gradient, and why the lab refuses to derive a contour direction by rotating a gradient. |
| 6 | Step size | Seven: the `h^2` law, the two error sources, and the best `h` for each of the two methods. |
| 7 | The zero gradient | Eight: three identical zero gradients, what they cannot distinguish, and what object you would need instead. |
| 8 | Models and cost | Nine: a three-parameter loss and its gradient by hand, the cost of a numerical gradient, and what autodiff changes. |

Work them out before running anything. A lab about derivatives whose answers
you can only get by running it is a lab that teaches you to trust output.

Two of these deserve a warning:

- **2.3** asks for `df/dx` of `xy` at the point `(1, 0)`. The answer is
  surprising, and 2.5 asks you what it does and does not imply.
- **4.6** asks whether a sweep of 360 directions will find *exactly* the
  largest possible rate of change. Think about what "360 directions" means
  before answering.

## What the numbers mean when you run it

```
206 passed                      finished
1 passed, 205 skipped           untouched
150 passed, 56 skipped          two thirds done, nothing wrong
149 passed, 1 failed, 56 skipped   one thing attempted and wrong
```

A failure is information, not a scolding. It prints your value and the correct
one side by side, and the assertion message usually names the specific mistake.

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

Almost every one of these can be done on paper. That is deliberate: a lab about
derivatives whose answers you cannot check by hand is a lab that teaches you to
trust output.

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

Check yourself from the LAB DIRECTORY:

    .venv/bin/pytest starter -q

Throughout, these are the six surfaces from `surfaces.py`:

    bowl     f = x^2 + 3y^2
    plane    f = 3x - 2y + 5
    product  f = xy
    saddle   f = x^2 - y^2
    dome     f = -(x^2 + y^2)
    cubic    f = x^3 + x*y^2
"""

# =============================================================================
# Exercise 2 -- partial derivatives by hand
# =============================================================================

# 2.1 f = x^2 + 3y^2. What is df/dx at the point (2, 1)? A float.
#     Freeze y at 1, differentiate x^2 + 3 with respect to x, put x = 2.
BOWL_DF_DX_AT_2_1 = None

# 2.2 The same function and the same point. What is df/dy? A float.
BOWL_DF_DY_AT_2_1 = None

# 2.3 f = xy. What is df/dx at (1, 0)? A float.
#     Careful: this is the one that catches people.
PRODUCT_DF_DX_AT_1_0 = None

# 2.4 f = xy at the same point (1, 0). What is df/dy? A float.
PRODUCT_DF_DY_AT_1_0 = None

# 2.5 Given your answers to 2.3 and 2.4: is the surface flat at (1, 0)?
#     Answer "yes" or "no".
IS_THE_PRODUCT_FLAT_AT_1_0 = None

# 2.6 f = x^3 + x*y^2. What is df/dx at (2, 1)? An integer or a float.
#     df/dx = 3x^2 + y^2. Substitute.
CUBIC_DF_DX_AT_2_1 = None

# 2.7 Same function, same point. What is df/dy? df/dy = 2xy.
CUBIC_DF_DY_AT_2_1 = None

# 2.8 Why is the symbol written with a rounded d rather than a straight one?
#     Answer with one of these strings:
#       "it is a different kind of derivative with different rules"
#       "it signals that the function has other inputs being held fixed"
#       "it means the answer is approximate rather than exact"
WHY_THE_ROUNDED_D = None


# =============================================================================
# Exercise 3 -- the gradient as a vector
# =============================================================================

# 3.1 The gradient of the bowl at (2, 1), as a list of two floats.
#     This is just 2.1 and 2.2 side by side.
BOWL_GRADIENT_AT_2_1 = None

# 3.2 How many components does the gradient of f(x, y) = x^2 + 3y^2 have?
#     An integer. The graph of this function is a surface in three dimensions;
#     the question is about the gradient, not the graph.
BOWL_GRADIENT_LENGTH = None

# 3.3 The gradient of the plane 3x - 2y + 5 at (100, -100), as a list of two
#     floats. Think before you compute anything.
PLANE_GRADIENT_FAR_AWAY = None

# 3.4 |grad f| for the bowl at (1, 1), where the gradient is (2, 6).
#     A float. This is Day 99's norm: sqrt(2^2 + 6^2).
BOWL_GRADIENT_MAGNITUDE_AT_1_1 = None

# 3.5 What does that magnitude MEAN, in words?
#     Answer with one of these strings:
#       "the height of the surface at that point"
#       "the rate of climb in the steepest direction, per unit of distance"
#       "the distance from the point to the minimum"
WHAT_THE_MAGNITUDE_MEANS = None

# 3.6 A model has 500 parameters. How many numbers are in the gradient of its
#     loss? An integer.
GRADIENT_LENGTH_FOR_500_PARAMETERS = None


# =============================================================================
# Exercise 4 -- directional derivatives and steepest ascent
# =============================================================================

# 4.1 At (1, 1) on the bowl the gradient is (2, 6). What is the directional
#     derivative along the direction (1, 0)? A float.
#     Remember that a directional derivative uses the UNIT direction, and
#     (1, 0) already has length 1.
BOWL_RATE_DUE_EAST = None

# 4.2 The same point, along the direction (0, 1). A float.
BOWL_RATE_DUE_NORTH = None

# 4.3 The same point, along the direction (3, -1). A float.
#     Work out the dot product of (2, 6) with (3, -1) first, then think about
#     what dividing by the length of (3, -1) does to a zero.
BOWL_RATE_ALONG_3_MINUS_1 = None

# 4.4 What is the LARGEST directional derivative available at that point, over
#     every possible unit direction? A float, to three decimal places or
#     better. You already computed it in 3.4.
BOWL_LARGEST_POSSIBLE_RATE = None

# 4.5 And the smallest (that is, the most negative)? A float.
BOWL_SMALLEST_POSSIBLE_RATE = None

# 4.6 If you sweep 360 directions one degree apart and take the largest rate,
#     will it exactly equal your answer to 4.4?
#     Answer with one of these strings:
#       "yes, exactly"
#       "no, slightly smaller"
#       "no, slightly larger"
WILL_THE_SWEEP_HIT_THE_MAXIMUM = None

# 4.7 A direction u makes an angle A with the gradient. The directional
#     derivative along u equals |grad f| times WHAT function of A?
#     Answer with one of these strings: "sin", "cos", "tan"
WHICH_TRIG_FUNCTION = None


# =============================================================================
# Exercise 5 -- contours and perpendicularity
# =============================================================================

# 5.1 A contour (or level set) of f is the set of points where f takes one
#     fixed value. On the bowl x^2 + 3y^2, what shape are the contours?
#     Answer with one of these strings: "circles", "ellipses", "straight lines"
BOWL_CONTOUR_SHAPE = None

# 5.2 On the plane 3x - 2y + 5, what shape are the contours?
#     Same three choices.
PLANE_CONTOUR_SHAPE = None

# 5.3 What is the angle, in degrees, between the gradient at a point and the
#     tangent to the contour through that point? A number.
ANGLE_BETWEEN_GRADIENT_AND_CONTOUR = None

# 5.4 If you walk a very short distance ALONG a contour, roughly how much does
#     f change?
#     Answer with one of these strings:
#       "it grows at the rate |grad f|"
#       "essentially nothing, to first order"
#       "it shrinks at the rate |grad f|"
WHAT_HAPPENS_ALONG_A_CONTOUR = None

# 5.5 The lab checks perpendicularity by taking two points a distance delta
#     apart on an exactly parametrised contour and dotting the chord between
#     them with the unit gradient. The answer is not exactly zero. When delta
#     is divided by 10, what happens to the dot product?
#     Answer with one of these strings:
#       "it stays the same"
#       "it is divided by about 10"
#       "it is divided by about 100"
HOW_THE_DOT_PRODUCT_SHRINKS = None

# 5.6 Why does the lab parametrise each contour algebraically instead of
#     finding the contour direction by rotating the gradient 90 degrees?
#     Answer with one of these strings:
#       "rotating is slower to compute"
#       "rotating would make the result true by construction and prove nothing"
#       "rotating only works in two dimensions"
WHY_NOT_ROTATE_THE_GRADIENT = None


# =============================================================================
# Exercise 6 -- step size, and Day 108's U-curve
# =============================================================================

# 6.1 For f = x^2, the central difference ((x+h)^2 - (x-h)^2) / (2h) simplifies
#     to what? Answer with one of these strings: "2x", "2x + h", "2x + h^2"
CENTRAL_DIFFERENCE_ON_A_SQUARE = None

# 6.2 For f = x^3, the same expression simplifies to 3x^2 plus what?
#     Answer with one of these strings: "0", "h", "h^2", "h^3"
CENTRAL_DIFFERENCE_ERROR_ON_A_CUBE = None

# 6.3 On a cubic, if you divide h by 10, the METHOD error is divided by what?
#     An integer.
TRUNCATION_ERROR_IMPROVEMENT_PER_DECADE = None

# 6.4 As h gets very small, a second source of error takes over. What is it?
#     Answer with one of these strings:
#       "the function becomes non-differentiable"
#       "subtracting two nearly equal floats loses the digits they shared"
#       "numpy switches to a lower precision"
WHAT_GOES_WRONG_FOR_TINY_H = None

# 6.5 Which h in the range 1e-1 down to 1e-14 gives the SMALLEST total error
#     for a central difference on the cubic, on float64? A float, such as
#     1e-05. The trough sits near the cube root of machine epsilon.
BEST_H_FOR_CENTRAL = None

# 6.6 And for a FORWARD difference, whose method error shrinks like h rather
#     than h^2, so the trough sits near the square root of machine epsilon?
#     A float.
BEST_H_FOR_FORWARD = None

# 6.7 At h = 1e-14 the central difference on the cubic gives an answer that is
#     compared with the answer at h = 0.1:
#     Answer with one of these strings: "much better", "about the same",
#     "much worse"
CENTRAL_AT_TINY_H_VERSUS_MODERATE_H = None


# =============================================================================
# Exercise 7 -- the zero gradient, and what it does not tell you
# =============================================================================

# 7.1 The gradient of the bowl at the origin, as a list of two floats.
BOWL_GRADIENT_AT_ORIGIN = None

# 7.2 The gradient of the saddle x^2 - y^2 at the origin, as a list of two
#     floats.
SADDLE_GRADIENT_AT_ORIGIN = None

# 7.3 The gradient of the dome -(x^2 + y^2) at the origin, as a list of two
#     floats.
DOME_GRADIENT_AT_ORIGIN = None

# 7.4 Given three identical answers above: can the gradient alone tell you
#     which of the three points is a minimum? Answer "yes" or "no".
CAN_THE_GRADIENT_TELL_THEM_APART = None

# 7.5 What is the general name for a point where every partial derivative is
#     zero? Answer with one of these strings:
#       "minimum", "stationary point", "inflection point"
NAME_FOR_A_ZERO_GRADIENT_POINT = None

# 7.6 Which object would you need in order to tell a minimum from a maximum
#     from a saddle -- the matrix of SECOND partial derivatives?
#     Answer with one of these strings: "Jacobian", "Hessian", "Laplacian"
WHAT_YOU_NEED_INSTEAD = None

# 7.7 On the saddle x^2 - y^2, walking 0.5 due east from the origin changes f
#     by how much? A float, with its sign.
SADDLE_CHANGE_WALKING_EAST = None

# 7.8 And 0.5 due north? A float, with its sign.
SADDLE_CHANGE_WALKING_NORTH = None


# =============================================================================
# Exercise 8 -- models, cost, and the AI thread
# =============================================================================

# 8.1 The loss L = (1/4) sum (w1*a + w2*b + c - y)^2 over the four samples in
#     surfaces.py, evaluated at w1 = w2 = c = 1. A float.
#     The four predictions are 4, 4, 7, 2 against targets 8, 7, 15, 3.
MODEL_LOSS_AT_ONES = None

# 8.2 dL/dw1 at that point. An integer or a float.
#     dL/dw1 = (2/4) sum (residual * a), and the four residuals are
#     -4, -3, -8, -1 with a values 1, 2, 3, 0.
MODEL_DL_DW1 = None

# 8.3 dL/dw2 there. The b values are 2, 1, 3, 1.
MODEL_DL_DW2 = None

# 8.4 dL/dc there. The c term has a coefficient of 1 in every sample.
MODEL_DL_DC = None

# 8.5 How many separate evaluations of the loss does ONE numerical gradient of
#     a 3-parameter model cost, using a central difference? An integer.
EVALUATIONS_FOR_A_3_PARAMETER_GRADIENT = None

# 8.6 And for a model with 1,000,000 parameters? An integer.
EVALUATIONS_FOR_A_MILLION_PARAMETER_GRADIENT = None

# 8.7 Reverse-mode automatic differentiation gets the whole gradient for a cost
#     that does what as the parameter count grows?
#     Answer with one of these strings:
#       "grows in proportion to the number of parameters"
#       "stays roughly one forward pass plus one backward pass"
#       "grows as the square of the number of parameters"
COST_OF_REVERSE_MODE_AUTODIFF = None

# 8.8 Numerical differentiation is still genuinely useful in training code, for
#     one specific job. Which?
#     Answer with one of these strings:
#       "it is faster than autodiff for small models"
#       "checking that a hand-written backward pass is correct"
#       "it handles non-differentiable functions that autodiff cannot"
WHAT_NUMERICAL_GRADIENTS_ARE_STILL_FOR = None

# 8.9 To make a loss go DOWN, you step along which vector?
#     Answer with one of these strings:
#       "the gradient"
#       "the negative gradient"
#       "any direction perpendicular to the gradient"
WHICH_WAY_TO_STEP_TO_REDUCE_A_LOSS = None

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

Both `examples/` and `starter/` contain modules called `gradients` and
`surfaces`, 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 `gradients` 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
`gradients`, `surfaces` 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 ("gradients", "surfaces", "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 -- eight functions to write. Your work goes here.

Each one raises NotImplementedError until you write it, and the test suite
SKIPS anything still unwritten rather than failing it. Your score only ever
counts work you actually attempted.

Check yourself from the LAB DIRECTORY (the one above this file):

    .venv/bin/pytest starter -q

Read each docstring before writing the body. Every one gives the formula and a
worked example small enough to check on paper.

The order matters. `partial` is the only one that touches f directly; every
other function here is built out of it, so if `partial` is wrong, everything
below inherits the error. Get it passing first.

The three helpers at the bottom are written for you.
"""

from __future__ import annotations

import numpy as np

from surfaces import H_DEFAULT, N_DIRECTIONS


# ===========================================================================
# The two that matter
# ===========================================================================


def partial(f, point, index, h=H_DEFAULT):
    """1.1 -- the partial derivative of f with respect to input `index`.

    Move coordinate `index` up by h and down by h, leave EVERY other
    coordinate exactly where it was, and divide the change in f by 2h:

        ( f(... x_i + h ...) - f(... x_i - h ...) ) / (2h)

    Three things to get right, in order of how often they are got wrong:

      * Copy the point before you modify it. `np.asarray(point, dtype=float)`
        then `.copy()` twice. If you mutate the caller's array, the tests will
        catch you, but your own code will be haunted.
      * Only ONE coordinate moves. That is the entire definition. A test
        watches which points f is actually called with.
      * Divide by 2h, not by h. Dividing by h gives an answer exactly twice
        too big, which is a satisfying bug because everything still looks
        plausible.

    Return a plain float, not a NumPy scalar: wrap the result in `float(...)`.

    >>> partial(lambda p: p[0] ** 2 + 3 * p[1] ** 2, (2.0, 1.0), 0)
    4.000000000026205
    """
    raise NotImplementedError("partial")


def gradient(f, point, h=H_DEFAULT):
    """1.2 -- the gradient: one partial derivative per input, as a vector.

    Call `partial` once for each coordinate of the point and put the results
    into a NumPy array, in order. Two or three lines.

    Return an array whose length matches the INPUT, not the number of
    dimensions the surface lives in: a function of two inputs has a
    two-component gradient, even though its graph is a surface in three
    dimensions.

    Use `np.asarray(point, dtype=float).size` to find out how many inputs
    there are rather than assuming two -- exercise 1.2 is tested on a function
    of three.

    >>> gradient(lambda p: p[0] ** 2 + 3 * p[1] ** 2, (1.0, 1.0)).round(6).tolist()
    [2.0, 6.0]
    """
    raise NotImplementedError("gradient")


# ===========================================================================
# Vector arithmetic -- Day 99, unchanged
# ===========================================================================


def magnitude(vector):
    """1.3 -- the length of a vector. Return a plain float.

    Day 99's Euclidean norm: the square root of the sum of the squares.
    `np.dot(v, v)` gives you that sum in one call.

    Applied to a gradient, this answers "how steep is the steepest way up".

    >>> magnitude([3.0, 4.0])
    5.0
    """
    raise NotImplementedError("magnitude")


def unit(vector):
    """1.4 -- the same direction, scaled to length exactly 1.

    Divide the vector by its magnitude. RAISE `ValueError` if the magnitude is
    zero, with a message containing the words "no direction" -- a zero vector
    has no direction to preserve, and returning NaNs instead of saying so
    would push the failure somewhere harder to find.

    >>> unit([3.0, 4.0]).tolist()
    [0.6, 0.8]
    """
    raise NotImplementedError("unit")


# ===========================================================================
# Directional derivatives -- Day 103's dot product, doing real work
# ===========================================================================


def directional_derivative(f, point, direction, h=H_DEFAULT):
    """1.5 -- how fast f changes if you walk from `point` along `direction`.

    Two lines. Normalise the direction with `unit`, then dot it with the
    gradient. Return a plain float.

    Normalising first is not tidiness. Without it, handing in an arrow twice
    as long would double the answer, and "the rate of change in this
    direction" would depend on how long an arrow you happened to draw.

    >>> f = lambda p: p[0] ** 2 + 3 * p[1] ** 2
    >>> round(directional_derivative(f, (1.0, 1.0), (1.0, 0.0)), 6)
    2.0
    """
    raise NotImplementedError("directional_derivative")


def directional_derivative_direct(f, point, direction, h=H_DEFAULT):
    """1.6 -- the same quantity, measured WITHOUT forming a gradient.

    Step h forward along the unit direction and h back along it, and divide by
    2h -- exactly `partial`, except that the step is along an arbitrary
    bearing instead of along an axis:

        ( f(p + h*u) - f(p - h*u) ) / (2h)

    No partials, no dot product, no assumption that the two routes agree.
    This function exists so that 1.5 can be CHECKED rather than believed, and
    the test that compares them is the most important one in the lab.

    >>> f = lambda p: p[0] ** 2 + 3 * p[1] ** 2
    >>> round(directional_derivative_direct(f, (1.0, 1.0), (0.0, 1.0)), 6)
    6.0
    """
    raise NotImplementedError("directional_derivative_direct")


def sweep_directions(f, point, n=None, h=H_DEFAULT):
    """1.7 -- try n bearings evenly spaced around the circle; report each rate.

    Return `(angles, rates)`, both NumPy arrays of length n:

      * `angles` -- n values from 0 up to but NOT including 2*pi. That is
        exactly `np.linspace(0.0, 2.0 * np.pi, n, endpoint=False)`; including
        the endpoint would sample 0 and 2*pi as two different bearings when
        they are the same one.
      * `rates` -- for each angle a, the rate of change along the direction
        `(cos a, sin a)`, measured with `directional_derivative_direct` so
        that the gradient plays no part in producing the numbers.

    Default n to N_DIRECTIONS, imported at the top of this file.

    Two-input functions only; this is the case you can draw.

    >>> f = lambda p: p[0] ** 2 + 3 * p[1] ** 2
    >>> angles, rates = sweep_directions(f, (1.0, 1.0), n=4)
    >>> rates.round(6).tolist()
    [2.0, 6.0, -2.0, -6.0]
    """
    raise NotImplementedError("sweep_directions")


def forward_partial(f, point, index, h=H_DEFAULT):
    """1.8 -- the one-sided version, for the comparison in exercise 6.

        ( f(x + h) - f(x) ) / h

    Return a plain float.

    One evaluation cheaper than a central difference when you already have
    f(x), and markedly worse: its error shrinks like h where a central
    difference's shrinks like h squared. You will measure exactly how much
    worse.

    >>> round(forward_partial(lambda p: p[0] ** 2, (2.0,), 0, 0.1), 6)
    4.1
    """
    raise NotImplementedError("forward_partial")


# ===========================================================================
# Written for you. Read them -- the tests use them.
# ===========================================================================


def angle_degrees(vector):
    """The bearing of a 2-D vector from the positive x-axis, in [0, 360)."""
    v = np.asarray(vector, dtype=float)
    return float(np.degrees(np.arctan2(v[1], v[0])) % 360.0)


def angular_gap_degrees(a, b):
    """The smaller of the two ways round between two bearings, in degrees.

    Without the wrap-around, 359.6 and 0.1 would look 359.5 degrees apart
    instead of 0.5, and every steepest-ascent check would fail for a reason
    that has nothing to do with calculus.
    """
    raw = abs(a - b) % 360.0
    return float(min(raw, 360.0 - raw))


def contour_chord(f, contour, level, t, delta):
    """A unit vector along a contour of f, built WITHOUT using the gradient.

    Takes two points on the exact algebraic contour, at parameters t and
    t + delta, and returns the unit vector from the first to the second, the
    two points, and f at each of them -- so the caller can check the curve
    really did stay on one level rather than trust the algebra.

    Note what this does not do: it never rotates the gradient. Deriving the
    contour direction from the gradient and then observing that the two are
    perpendicular would prove nothing at all.
    """
    p = contour(level, t)
    q = contour(level, t + delta)
    return unit(q - p), p, q, f(p), f(q)

"""The surfaces this lab measures, and the exact gradients to check against.

Every function here is a function of several inputs. Every one has a gradient
you can work out with a pencil in under a minute, which is the entire point:
the numerical machinery in `gradients.py` is only trustworthy if there is
something exact to hold it against.

All the data is invented. None of it is a measurement of anything real. What
IS real is every number the lab prints about these functions, because those are
computed from these definitions at run time.

Read this file. Do not change it -- the tests compare against the values
written down here.
"""

from __future__ import annotations

import numpy as np

# --------------------------------------------------------------------------
# The step size, and why this one
# --------------------------------------------------------------------------
#
# Day 108 established the shape: a central difference has a truncation error
# that shrinks like h squared and a roundoff error that GROWS like 1/h, so the
# total error is U-shaped in h and the best step is somewhere in the middle.
# For a central difference on float64 the trough sits near the cube root of
# the machine epsilon, which is about 6e-6. This lab uses 1e-5, which is close
# enough to the bottom of that trough to be within a factor of two of the best
# achievable error on every surface here, and is a round number a reader can
# remember. Script 06 sweeps h over twelve orders of magnitude and prints the
# curve rather than asking you to take this on trust.
H_DEFAULT = 1.0e-5

# The tolerance every gradient assertion in this lab uses.
#
# It is set from what the arithmetic can achieve, not from what makes the
# tests pass. Four of the six surfaces below are at most quadratic in each
# variable, and a central difference is ALGEBRAICALLY EXACT for those -- the
# h-squared term is multiplied by a third derivative that is zero, so the only
# error left is floating-point roundoff, which lands around 1e-11 at h = 1e-5
# for points of this size. The one genuinely cubic surface has a truncation
# error of exactly h squared, which is 1e-10 here. GRADIENT_TOL is three
# orders of magnitude above the worst of those, which leaves room for a
# different processor's rounding without leaving room for a wrong answer.
GRADIENT_TOL = 1.0e-8

# Directions sampled around the full circle when the lab checks that the
# gradient really is the steepest way up. With 360 evenly spaced directions
# the closest sample sits within 0.5 degrees of any given angle, so the
# tolerance below is that bound with a little slack for the wrap-around at
# 360 degrees.
N_DIRECTIONS = 360
ANGLE_TOL_DEGREES = 1.0

# The step taken along a contour when the lab checks perpendicularity. The
# chord between two points on a curve is not the tangent; it differs from it
# by an angle of order delta, so the dot product of a unit gradient with a
# unit chord is of order delta rather than exactly zero. Script 04 halves
# delta four times and prints the dot product each time, so the reader watches
# it shrink instead of being handed a tolerance.
CONTOUR_DELTA = 1.0e-5
CONTOUR_DOT_TOL = 1.0e-4

# Used where a random point or a random direction is wanted. Seeded, so every
# number in `expected-output/` is reproducible.
SEED = 109


# --------------------------------------------------------------------------
# 1. A quadratic bowl -- the shape every optimisation picture is drawn on
# --------------------------------------------------------------------------

def bowl(point):
    """f(x, y) = x^2 + 3y^2. A bowl with its lowest point at the origin.

    The 3 makes it an elliptical bowl rather than a circular one: it is three
    times steeper in y than in x at the same distance out, which is exactly
    the situation that makes gradient descent zig-zag on Day 111.
    """
    x, y = point
    return x * x + 3.0 * y * y


def bowl_gradient(point):
    """grad f = (2x, 6y). Differentiate x^2 + 3y^2 one variable at a time."""
    x, y = point
    return np.array([2.0 * x, 6.0 * y])


# --------------------------------------------------------------------------
# 2. A plane -- the function whose gradient is the same everywhere
# --------------------------------------------------------------------------

def plane(point):
    """f(x, y) = 3x - 2y + 5. A flat tilted sheet."""
    x, y = point
    return 3.0 * x - 2.0 * y + 5.0


def plane_gradient(point):
    """grad f = (3, -2), at every point in the plane. The point is ignored."""
    del point
    return np.array([3.0, -2.0])


# --------------------------------------------------------------------------
# 3. A product -- the smallest function whose partials involve each other
# --------------------------------------------------------------------------

def product(point):
    """f(x, y) = xy. Flat along both axes through the origin, curved between.

    This is the function that makes the phrase "hold the other one still" do
    real work. Along the x-axis (y = 0) the function is identically zero, so
    the slope in x is zero. Move off that line and the slope in x is y.
    """
    x, y = point
    return x * y


def product_gradient(point):
    """grad f = (y, x). Each partial is the OTHER variable, held fixed."""
    x, y = point
    return np.array([y, x])


# --------------------------------------------------------------------------
# 4. A saddle -- a stationary point that is neither a peak nor a floor
# --------------------------------------------------------------------------

def saddle(point):
    """f(x, y) = x^2 - y^2. Up along x, down along y, flat at the origin."""
    x, y = point
    return x * x - y * y


def saddle_gradient(point):
    """grad f = (2x, -2y). Zero at the origin, and the origin is a saddle."""
    x, y = point
    return np.array([2.0 * x, -2.0 * y])


# --------------------------------------------------------------------------
# 5. A dome -- a stationary point that IS a maximum
# --------------------------------------------------------------------------

def dome(point):
    """f(x, y) = -(x^2 + y^2). The bowl turned upside down."""
    x, y = point
    return -(x * x + y * y)


def dome_gradient(point):
    """grad f = (-2x, -2y). Also zero at the origin. Also a stationary point."""
    x, y = point
    return np.array([-2.0 * x, -2.0 * y])


# --------------------------------------------------------------------------
# 6. A genuine cubic -- the one surface where truncation error is visible
# --------------------------------------------------------------------------

def cubic(point):
    """f(x, y) = x^3 + x*y^2.

    Every other surface here is at most quadratic in each variable, which
    makes the central difference exact for them. This one is not, and its
    error is not merely small but PREDICTABLE: expanding
    ((x+h)^3 - (x-h)^3) / (2h) gives 3x^2 + h^2 exactly, so the numerical
    partial in x overshoots the true one by exactly h squared, with no other
    terms at all. Script 06 measures that and checks it to twelve decimal
    places.
    """
    x, y = point
    return x ** 3 + x * y * y


def cubic_gradient(point):
    """grad f = (3x^2 + y^2, 2xy)."""
    x, y = point
    return np.array([3.0 * x * x + y * y, 2.0 * x * y])


# --------------------------------------------------------------------------
# The registry the scripts and tests iterate over
# --------------------------------------------------------------------------

SURFACES = {
    "bowl": (bowl, bowl_gradient, "x^2 + 3y^2", "grad = (2x, 6y)"),
    "plane": (plane, plane_gradient, "3x - 2y + 5", "grad = (3, -2)"),
    "product": (product, product_gradient, "xy", "grad = (y, x)"),
    "saddle": (saddle, saddle_gradient, "x^2 - y^2", "grad = (2x, -2y)"),
    "dome": (dome, dome_gradient, "-(x^2 + y^2)", "grad = (-2x, -2y)"),
    "cubic": (cubic, cubic_gradient, "x^3 + x*y^2", "grad = (3x^2 + y^2, 2xy)"),
}

# The points every surface is probed at. Chosen by hand: one in the first
# quadrant, one with a negative coordinate, one off-axis with a fraction, one
# ON an axis where a partial vanishes, and one far out.
PROBE_POINTS = (
    (1.0, 1.0),
    (2.0, -1.0),
    (-0.5, 3.0),
    (4.0, 0.0),
    (0.25, 0.75),
)

# The three surfaces that all have a zero gradient at the origin, and what the
# origin actually IS for each. The gradient cannot tell them apart; this table
# is the answer the gradient does not carry.
STATIONARY_AT_ORIGIN = (
    ("bowl", "minimum", "every direction goes up"),
    ("dome", "maximum", "every direction goes down"),
    ("saddle", "saddle", "up along x, down along y"),
)


# --------------------------------------------------------------------------
# Exact contours, parametrised WITHOUT reference to the gradient
# --------------------------------------------------------------------------
#
# The claim being tested is that the gradient is perpendicular to the contour.
# Walking along the contour by stepping perpendicular to the gradient would
# make that claim true by construction and prove nothing. So each contour
# below is an exact algebraic parametrisation, derived on paper from the
# function alone, and the lab checks it lands back on the same value of f
# before it uses it for anything.

def bowl_contour(level, t):
    """A point on the ellipse x^2 + 3y^2 = level, for the angle parameter t.

    Substitute x = sqrt(level) cos t and y = sqrt(level/3) sin t into
    x^2 + 3y^2 and you get level (cos^2 t + sin^2 t) = level, for every t.
    """
    a = np.sqrt(level)
    b = np.sqrt(level / 3.0)
    return np.array([a * np.cos(t), b * np.sin(t)])


def product_contour(level, t):
    """A point on the hyperbola xy = level, for the parameter t (x = t)."""
    return np.array([t, level / t])


def dome_contour(level, t):
    """A point on the circle -(x^2 + y^2) = level, for the angle t.

    level must be negative; the radius is sqrt(-level).
    """
    r = np.sqrt(-level)
    return np.array([r * np.cos(t), r * np.sin(t)])


CONTOURS = {
    "bowl": (bowl, bowl_contour, 4.0, 0.7),
    "product": (product, product_contour, 6.0, 2.0),
    "dome": (dome, dome_contour, -9.0, 1.1),
}


# --------------------------------------------------------------------------
# A three-parameter model, so "one partial per parameter" is not just a claim
# --------------------------------------------------------------------------
#
# Four invented samples. Each row is (a, b, target). Nothing was measured to
# produce these; they were chosen so the arithmetic stays checkable by hand.

SAMPLES = (
    (1.0, 2.0, 8.0),
    (2.0, 1.0, 7.0),
    (3.0, 3.0, 15.0),
    (0.0, 1.0, 3.0),
)

# The parameter vector the lab evaluates the loss and its gradient at.
START_PARAMS = (1.0, 1.0, 1.0)


def model_loss(params):
    """Mean squared error of `pred = w1*a + w2*b + c` over SAMPLES.

    Three parameters, so the gradient is a vector of three numbers -- one per
    parameter. A real network has millions of parameters and the gradient has
    millions of entries. The idea does not change; only the count does.
    """
    w1, w2, c = params
    total = 0.0
    for a, b, target in SAMPLES:
        residual = w1 * a + w2 * b + c - target
        total += residual * residual
    return total / len(SAMPLES)


def model_loss_gradient(params):
    """The exact gradient, differentiated by hand.

    L = (1/n) sum r_i^2 with r_i = w1*a_i + w2*b_i + c - y_i, so
    dL/dw1 = (2/n) sum r_i * a_i, dL/dw2 = (2/n) sum r_i * b_i,
    dL/dc  = (2/n) sum r_i.
    """
    w1, w2, c = params
    n = len(SAMPLES)
    g_w1 = g_w2 = g_c = 0.0
    for a, b, target in SAMPLES:
        residual = w1 * a + w2 * b + c - target
        g_w1 += 2.0 * residual * a
        g_w2 += 2.0 * residual * b
        g_c += 2.0 * residual
    return np.array([g_w1 / n, g_w2 / n, g_c / n])

"""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 message
prints both your answer and the real one.

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

import math

import numpy as np
import pytest

import answers
import surfaces as S
from gradients import (
    angle_degrees,
    angular_gap_degrees,
    contour_chord,
    directional_derivative,
    directional_derivative_direct,
    forward_partial,
    gradient,
    magnitude,
    partial,
    sweep_directions,
    unit,
)


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


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


def test_0_the_environment_is_ready():
    """Always passes once the install worked. Everything below is your work."""
    assert int(np.__version__.split(".")[0]) >= 2, "numpy 2 or later is importable"
    assert S.bowl((2.0, 1.0)) == 7.0, "surfaces.py loads and the bowl is the bowl"
    assert S.model_loss(S.START_PARAMS) == 22.5, "the model data is intact"


# -- Exercise 1: your gradients.py -------------------------------------------


def test_1_1_partial_on_the_bowl_in_x():
    got = written(lambda: partial(S.bowl, (2.0, 1.0), 0))
    assert got == pytest.approx(4.0, abs=S.GRADIENT_TOL)


def test_1_1_partial_on_the_bowl_in_y():
    got = written(lambda: partial(S.bowl, (2.0, 1.0), 1))
    assert got == pytest.approx(6.0, abs=S.GRADIENT_TOL)


def test_1_1_partial_divides_by_two_h_not_by_h():
    """The most common bug: an answer exactly twice too big."""
    got = written(lambda: partial(S.plane, (0.0, 0.0), 0))
    assert got != pytest.approx(6.0, abs=1e-6), "divide by 2h, not by h"
    assert got == pytest.approx(3.0, abs=S.GRADIENT_TOL)


def test_1_1_partial_holds_the_other_coordinate_completely_still():
    seen = []

    def spy(p):
        seen.append(tuple(float(v) for v in p))
        return S.product(p)

    written(lambda: partial(spy, (2.0, 5.0), 0))
    assert len(seen) == 2, "a central difference evaluates f exactly twice"
    assert {p[1] for p in seen} == {5.0}, "y must be identical in both calls"


def test_1_1_partial_does_not_mutate_the_point_it_was_given():
    point = np.array([1.0, 2.0])
    written(lambda: partial(S.bowl, point, 0))
    assert point.tolist() == [1.0, 2.0], "copy the point before nudging it"


def test_1_1_partial_returns_a_plain_float():
    got = written(lambda: partial(S.bowl, (1.0, 1.0), 0))
    assert isinstance(got, float), "wrap the result in float(...)"


@pytest.mark.parametrize("name", sorted(S.SURFACES))
@pytest.mark.parametrize("point", S.PROBE_POINTS)
@pytest.mark.parametrize("index", (0, 1))
def test_1_1_partial_on_every_surface_and_point(name, point, index):
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    got = written(lambda: partial(f, point, index))
    assert got == pytest.approx(exact_gradient(point)[index], abs=S.GRADIENT_TOL)


def test_1_2_gradient_of_the_bowl():
    got = written(lambda: gradient(S.bowl, (1.0, 1.0)))
    assert got == pytest.approx([2.0, 6.0], abs=S.GRADIENT_TOL)


def test_1_2_gradient_returns_a_numpy_array():
    got = written(lambda: gradient(S.bowl, (1.0, 1.0)))
    assert isinstance(got, np.ndarray), "return an array, not a list or a tuple"


def test_1_2_gradient_has_one_entry_per_INPUT():
    got = written(lambda: gradient(S.bowl, (1.0, 1.0)))
    assert got.size == 2, "two inputs, two partials -- not three"


def test_1_2_gradient_works_on_a_function_of_three_inputs():
    got = written(lambda: gradient(S.model_loss, S.START_PARAMS))
    assert got.size == 3, "do not hard-code two coordinates"
    assert got == pytest.approx([-17.0, -18.0, -8.0], abs=S.GRADIENT_TOL)


@pytest.mark.parametrize("name", sorted(S.SURFACES))
@pytest.mark.parametrize("point", S.PROBE_POINTS)
def test_1_2_gradient_on_every_surface_and_point(name, point):
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    got = written(lambda: gradient(f, point))
    assert got == pytest.approx(exact_gradient(point), abs=S.GRADIENT_TOL)


def test_1_3_magnitude_of_a_three_four_five_triangle():
    assert written(lambda: magnitude([3.0, 4.0])) == 5.0


def test_1_3_magnitude_of_the_zero_vector_is_zero():
    assert written(lambda: magnitude([0.0, 0.0])) == 0.0


def test_1_3_magnitude_returns_a_plain_float():
    assert isinstance(written(lambda: magnitude([1.0, 1.0])), float)


def test_1_3_magnitude_of_the_bowl_gradient_at_one_one():
    got = written(lambda: magnitude(gradient(S.bowl, (1.0, 1.0))))
    assert got == pytest.approx(math.sqrt(40.0), abs=1e-9)


def test_1_4_unit_scales_to_length_one():
    got = written(lambda: unit([3.0, 4.0]))
    assert got == pytest.approx([0.6, 0.8], abs=1e-15)


def test_1_4_unit_keeps_the_bearing():
    got = written(lambda: unit([2.0, 6.0]))
    assert angle_degrees(got) == pytest.approx(angle_degrees([2.0, 6.0]), abs=1e-9)


def test_1_4_unit_of_a_long_and_a_short_arrow_agree():
    a = written(lambda: unit([1.0, 2.0]))
    b = written(lambda: unit([50.0, 100.0]))
    assert a == pytest.approx(b, abs=1e-15)


def test_1_4_unit_raises_value_error_on_the_zero_vector():
    def call():
        with pytest.raises(ValueError, match="no direction"):
            unit([0.0, 0.0])
        return True

    assert written(call) is True


def test_1_5_directional_derivative_due_east_is_the_x_partial():
    got = written(lambda: directional_derivative(S.bowl, (1.0, 1.0), (1.0, 0.0)))
    assert got == pytest.approx(2.0, abs=S.GRADIENT_TOL)


def test_1_5_directional_derivative_due_north_is_the_y_partial():
    got = written(lambda: directional_derivative(S.bowl, (1.0, 1.0), (0.0, 1.0)))
    assert got == pytest.approx(6.0, abs=S.GRADIENT_TOL)


def test_1_5_directional_derivative_normalises_the_direction():
    """A longer arrow must not give a bigger answer."""
    short = written(lambda: directional_derivative(S.bowl, (1.0, 1.0), (1.0, 2.0)))
    long = written(lambda: directional_derivative(S.bowl, (1.0, 1.0), (1000.0, 2000.0)))
    assert short == pytest.approx(long, abs=1e-9), "normalise before dotting"


def test_1_5_a_direction_perpendicular_to_the_gradient_gives_zero():
    got = written(lambda: directional_derivative(S.bowl, (1.0, 1.0), (3.0, -1.0)))
    assert got == pytest.approx(0.0, abs=S.GRADIENT_TOL)


def test_1_6_direct_measurement_agrees_with_the_dot_product():
    """The most important test in the lab: two routes, one answer."""
    for direction in ((1.0, 1.0), (-1.0, 2.0), (7.0, 0.5), (-2.0, -5.0)):
        via = written(lambda d=direction: directional_derivative(S.bowl, (1.0, 1.0), d))
        direct = written(
            lambda d=direction: directional_derivative_direct(S.bowl, (1.0, 1.0), d)
        )
        assert via == pytest.approx(direct, abs=S.GRADIENT_TOL)


def test_1_6_direct_measurement_never_forms_a_gradient():
    """f must be evaluated exactly twice, not once per axis plus twice more."""
    calls = []

    def counted(p):
        calls.append(1)
        return S.bowl(p)

    written(lambda: directional_derivative_direct(counted, (1.0, 1.0), (1.0, 1.0)))
    assert len(calls) == 2, "two evaluations: one forward, one back"


def test_1_6_walking_backwards_negates_the_rate():
    forward = written(
        lambda: directional_derivative_direct(S.cubic, (2.0, -1.0), (1.0, 3.0))
    )
    backward = written(
        lambda: directional_derivative_direct(S.cubic, (2.0, -1.0), (-1.0, -3.0))
    )
    assert forward == pytest.approx(-backward, abs=1e-9)


def test_1_7_sweep_returns_two_arrays_of_the_right_length():
    angles, rates = written(lambda: sweep_directions(S.bowl, (1.0, 1.0), n=8))
    assert len(angles) == 8 and len(rates) == 8


def test_1_7_sweep_excludes_the_endpoint():
    """Sampling both 0 and 2*pi would count the same bearing twice."""
    angles, _rates = written(lambda: sweep_directions(S.bowl, (1.0, 1.0), n=4))
    assert float(angles[0]) == 0.0
    assert float(angles[-1]) == pytest.approx(1.5 * math.pi, abs=1e-12)


def test_1_7_sweep_of_four_gives_the_two_partials_and_their_negatives():
    _angles, rates = written(lambda: sweep_directions(S.bowl, (1.0, 1.0), n=4))
    assert rates == pytest.approx([2.0, 6.0, -2.0, -6.0], abs=S.GRADIENT_TOL)


def test_1_7_sweep_defaults_to_the_documented_number_of_directions():
    _angles, rates = written(lambda: sweep_directions(S.bowl, (1.0, 1.0)))
    assert len(rates) == S.N_DIRECTIONS


def test_1_8_forward_partial_on_a_square():
    got = written(lambda: forward_partial(S.bowl, (2.0, 1.0), 0, 0.1))
    assert got == pytest.approx(4.1, abs=1e-9), "((2.1)^2 - 2^2) / 0.1 = 4.1"


def test_1_8_forward_partial_evaluates_f_only_twice():
    calls = []

    def counted(p):
        calls.append(1)
        return S.bowl(p)

    written(lambda: forward_partial(counted, (1.0, 1.0), 0))
    assert len(calls) == 2


def test_1_8_forward_is_worse_than_central_at_the_default_step():
    exact = float(S.cubic_gradient((2.0, 1.0))[0])
    c = written(lambda: abs(partial(S.cubic, (2.0, 1.0), 0) - exact))
    f = written(lambda: abs(forward_partial(S.cubic, (2.0, 1.0), 0) - exact))
    assert f > 1000 * c


# -- Exercise 1 (applied): the two facts the day exists for -------------------


@pytest.mark.parametrize("name,point", (
    ("bowl", (1.0, 1.0)),
    ("bowl", (0.25, 0.75)),
    ("product", (2.0, -1.0)),
    ("saddle", (1.5, 0.5)),
    ("cubic", (1.0, 1.0)),
))
def test_applied_the_gradient_wins_the_sweep(name, point):
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    angles, rates = written(lambda: sweep_directions(f, point))
    best = int(np.argmax(rates))
    gap = angular_gap_degrees(
        float(np.degrees(angles[best])), angle_degrees(exact_gradient(point))
    )
    assert gap <= S.ANGLE_TOL_DEGREES


@pytest.mark.parametrize("name,point", (
    ("bowl", (1.0, 1.0)),
    ("product", (2.0, -1.0)),
    ("cubic", (1.0, 1.0)),
))
def test_applied_no_direction_beats_the_gradients_magnitude(name, point):
    f, exact_gradient = S.SURFACES[name][0], S.SURFACES[name][1]
    _angles, rates = written(lambda: sweep_directions(f, point))
    assert float(np.max(rates)) <= magnitude(exact_gradient(point)) + S.GRADIENT_TOL


@pytest.mark.parametrize("name", sorted(S.CONTOURS))
@pytest.mark.parametrize("t", (0.4, 0.9, 1.4, 1.9))
def test_applied_the_gradient_is_perpendicular_to_the_contour(name, t):
    f, contour, level, _t0 = S.CONTOURS[name]
    result = written(lambda: contour_chord(f, contour, level, t, S.CONTOUR_DELTA))
    chord, p = result[0], result[1]
    g = written(lambda: unit(gradient(f, p)))
    assert abs(float(np.dot(g, chord))) < S.CONTOUR_DOT_TOL


@pytest.mark.parametrize("name", sorted(S.CONTOURS))
def test_applied_the_dot_product_shrinks_with_the_step(name):
    f, contour, level, t0 = S.CONTOURS[name]
    previous = None
    for k in (2, 3, 4, 5):
        delta = 10.0 ** (-k)
        result = written(lambda d=delta: contour_chord(f, contour, level, t0, d))
        chord, p = result[0], result[1]
        dot = abs(float(np.dot(written(lambda: unit(gradient(f, p))), chord)))
        if previous is not None:
            assert 9.0 < previous / dot < 11.0
        previous = dot


def test_applied_the_plane_has_the_same_gradient_everywhere():
    seen = [written(lambda p=p: gradient(S.plane, p))
            for p in ((0.0, 0.0), (1.0, 1.0), (-40.0, 17.5))]
    assert float(np.max(np.abs(np.array(seen) - seen[0]))) < S.GRADIENT_TOL


def test_applied_the_bowls_gradient_points_away_from_its_minimum():
    for point in ((0.5, 0.5), (1.0, 1.0), (2.0, 2.0), (-3.0, 1.0)):
        g = written(lambda p=point: unit(gradient(S.bowl, p)))
        assert float(np.dot(g, unit(np.array(point)))) > 0.0


@pytest.mark.parametrize("name,_kind,_why", S.STATIONARY_AT_ORIGIN)
def test_applied_all_three_stationary_points_have_a_zero_gradient(name, _kind, _why):
    g = written(lambda: gradient(S.SURFACES[name][0], (0.0, 0.0)))
    assert magnitude(g) < S.GRADIENT_TOL


def test_applied_the_cubics_truncation_error_is_exactly_h_squared():
    exact = float(S.cubic_gradient((2.0, 1.0))[0])
    for h in (1e-1, 1e-2, 1e-3):
        error = written(lambda hh=h: partial(S.cubic, (2.0, 1.0), 0, hh) - exact)
        assert error == pytest.approx(h * h, rel=1e-5)


def test_applied_the_error_curve_is_u_shaped():
    exact = float(S.cubic_gradient((2.0, 1.0))[0])
    errors = {k: written(lambda kk=k: abs(partial(S.cubic, (2.0, 1.0), 0, 10.0 ** -kk) - exact))
              for k in range(0, 15)}
    best = min(errors, key=errors.get)
    assert 0 < best < 14
    assert 10.0 ** -best == S.H_DEFAULT


def test_applied_a_numerical_gradient_costs_two_evaluations_per_parameter():
    calls = []

    def counted(p):
        calls.append(1)
        return S.model_loss(p)

    written(lambda: gradient(counted, S.START_PARAMS))
    assert len(calls) == 6


def test_applied_a_step_against_the_gradient_reduces_the_loss():
    g = written(lambda: gradient(S.model_loss, S.START_PARAMS))
    before = S.model_loss(S.START_PARAMS)
    assert S.model_loss(np.array(S.START_PARAMS) - 0.01 * g) < before


# -- Exercise 2: partial derivatives by hand ---------------------------------


def test_2_1_bowl_df_dx():
    assert predicted("BOWL_DF_DX_AT_2_1") == pytest.approx(4.0)


def test_2_2_bowl_df_dy():
    assert predicted("BOWL_DF_DY_AT_2_1") == pytest.approx(6.0)


def test_2_3_product_df_dx_at_1_0():
    assert predicted("PRODUCT_DF_DX_AT_1_0") == pytest.approx(0.0)


def test_2_4_product_df_dy_at_1_0():
    assert predicted("PRODUCT_DF_DY_AT_1_0") == pytest.approx(1.0)


def test_2_5_the_product_is_not_flat_there():
    assert predicted("IS_THE_PRODUCT_FLAT_AT_1_0") == "no"


def test_2_6_cubic_df_dx():
    assert predicted("CUBIC_DF_DX_AT_2_1") == pytest.approx(13.0)


def test_2_7_cubic_df_dy():
    assert predicted("CUBIC_DF_DY_AT_2_1") == pytest.approx(4.0)


def test_2_8_why_the_rounded_d():
    assert predicted("WHY_THE_ROUNDED_D") == (
        "it signals that the function has other inputs being held fixed"
    )


# -- Exercise 3: the gradient as a vector ------------------------------------


def test_3_1_bowl_gradient():
    assert predicted("BOWL_GRADIENT_AT_2_1") == pytest.approx([4.0, 6.0])


def test_3_2_gradient_length():
    assert predicted("BOWL_GRADIENT_LENGTH") == 2


def test_3_3_plane_gradient_far_away():
    assert predicted("PLANE_GRADIENT_FAR_AWAY") == pytest.approx([3.0, -2.0])


def test_3_4_bowl_gradient_magnitude():
    assert predicted("BOWL_GRADIENT_MAGNITUDE_AT_1_1") == pytest.approx(
        math.sqrt(40.0), abs=1e-3
    )


def test_3_5_what_the_magnitude_means():
    assert predicted("WHAT_THE_MAGNITUDE_MEANS") == (
        "the rate of climb in the steepest direction, per unit of distance"
    )


def test_3_6_one_partial_per_parameter():
    assert predicted("GRADIENT_LENGTH_FOR_500_PARAMETERS") == 500


# -- Exercise 4: directional derivatives -------------------------------------


def test_4_1_rate_due_east():
    assert predicted("BOWL_RATE_DUE_EAST") == pytest.approx(2.0)


def test_4_2_rate_due_north():
    assert predicted("BOWL_RATE_DUE_NORTH") == pytest.approx(6.0)


def test_4_3_rate_along_a_perpendicular():
    assert predicted("BOWL_RATE_ALONG_3_MINUS_1") == pytest.approx(0.0, abs=1e-9)


def test_4_4_largest_possible_rate():
    assert predicted("BOWL_LARGEST_POSSIBLE_RATE") == pytest.approx(
        math.sqrt(40.0), abs=1e-3
    )


def test_4_5_smallest_possible_rate():
    assert predicted("BOWL_SMALLEST_POSSIBLE_RATE") == pytest.approx(
        -math.sqrt(40.0), abs=1e-3
    )


def test_4_6_the_sweep_falls_just_short():
    assert predicted("WILL_THE_SWEEP_HIT_THE_MAXIMUM") == "no, slightly smaller"


def test_4_7_which_trig_function():
    assert predicted("WHICH_TRIG_FUNCTION") == "cos"


# -- Exercise 5: contours ----------------------------------------------------


def test_5_1_bowl_contours_are_ellipses():
    assert predicted("BOWL_CONTOUR_SHAPE") == "ellipses"


def test_5_2_plane_contours_are_straight_lines():
    assert predicted("PLANE_CONTOUR_SHAPE") == "straight lines"


def test_5_3_the_angle_is_ninety_degrees():
    assert predicted("ANGLE_BETWEEN_GRADIENT_AND_CONTOUR") == pytest.approx(90.0)


def test_5_4_nothing_happens_along_a_contour():
    assert predicted("WHAT_HAPPENS_ALONG_A_CONTOUR") == (
        "essentially nothing, to first order"
    )


def test_5_5_the_dot_product_is_first_order_in_delta():
    assert predicted("HOW_THE_DOT_PRODUCT_SHRINKS") == "it is divided by about 10"


def test_5_6_why_not_rotate_the_gradient():
    assert predicted("WHY_NOT_ROTATE_THE_GRADIENT") == (
        "rotating would make the result true by construction and prove nothing"
    )


# -- Exercise 6: step size ---------------------------------------------------


def test_6_1_central_difference_on_a_square():
    assert predicted("CENTRAL_DIFFERENCE_ON_A_SQUARE") == "2x"


def test_6_2_central_difference_error_on_a_cube():
    assert predicted("CENTRAL_DIFFERENCE_ERROR_ON_A_CUBE") == "h^2"


def test_6_3_second_order_means_a_hundredfold():
    assert predicted("TRUNCATION_ERROR_IMPROVEMENT_PER_DECADE") == 100


def test_6_4_what_goes_wrong_for_tiny_h():
    assert predicted("WHAT_GOES_WRONG_FOR_TINY_H") == (
        "subtracting two nearly equal floats loses the digits they shared"
    )


def test_6_5_best_h_for_central():
    assert predicted("BEST_H_FOR_CENTRAL") == pytest.approx(1e-05)


def test_6_6_best_h_for_forward():
    assert predicted("BEST_H_FOR_FORWARD") == pytest.approx(1e-08)


def test_6_7_tiny_h_is_much_worse():
    assert predicted("CENTRAL_AT_TINY_H_VERSUS_MODERATE_H") == "much worse"


# -- Exercise 7: the zero gradient -------------------------------------------


def test_7_1_bowl_gradient_at_origin():
    assert predicted("BOWL_GRADIENT_AT_ORIGIN") == pytest.approx([0.0, 0.0])


def test_7_2_saddle_gradient_at_origin():
    assert predicted("SADDLE_GRADIENT_AT_ORIGIN") == pytest.approx([0.0, 0.0])


def test_7_3_dome_gradient_at_origin():
    assert predicted("DOME_GRADIENT_AT_ORIGIN") == pytest.approx([0.0, 0.0])


def test_7_4_the_gradient_cannot_tell_them_apart():
    assert predicted("CAN_THE_GRADIENT_TELL_THEM_APART") == "no"


def test_7_5_the_name_is_stationary_point():
    assert predicted("NAME_FOR_A_ZERO_GRADIENT_POINT") == "stationary point"


def test_7_6_you_would_need_the_hessian():
    assert predicted("WHAT_YOU_NEED_INSTEAD") == "Hessian"


def test_7_7_the_saddle_rises_going_east():
    assert predicted("SADDLE_CHANGE_WALKING_EAST") == pytest.approx(0.25)


def test_7_8_and_falls_going_north():
    assert predicted("SADDLE_CHANGE_WALKING_NORTH") == pytest.approx(-0.25)


# -- Exercise 8: models and cost ---------------------------------------------


def test_8_1_model_loss():
    assert predicted("MODEL_LOSS_AT_ONES") == pytest.approx(22.5)


def test_8_2_dl_dw1():
    assert predicted("MODEL_DL_DW1") == pytest.approx(-17.0)


def test_8_3_dl_dw2():
    assert predicted("MODEL_DL_DW2") == pytest.approx(-18.0)


def test_8_4_dl_dc():
    assert predicted("MODEL_DL_DC") == pytest.approx(-8.0)


def test_8_5_six_evaluations():
    assert predicted("EVALUATIONS_FOR_A_3_PARAMETER_GRADIENT") == 6


def test_8_6_two_million_evaluations():
    assert predicted("EVALUATIONS_FOR_A_MILLION_PARAMETER_GRADIENT") == 2_000_000


def test_8_7_reverse_mode_does_not_scale_with_parameters():
    assert predicted("COST_OF_REVERSE_MODE_AUTODIFF") == (
        "stays roughly one forward pass plus one backward pass"
    )


def test_8_8_gradient_checking():
    assert predicted("WHAT_NUMERICAL_GRADIENTS_ARE_STILL_FOR") == (
        "checking that a hand-written backward pass is correct"
    )


def test_8_9_step_against_the_gradient():
    assert predicted("WHICH_WAY_TO_STEP_TO_REDUCE_A_LOSS") == "the negative gradient"

#!/usr/bin/env bash
# Tests for the Day 109 lab. Run from the lab directory:
#   bash tests/run_tests.sh
#
# The harness proves the lesson's claims by running code and reading real
# values, never by reading source:
#
#   * a partial derivative moves ONE coordinate and holds the rest still, and
#     the harness watches which points f is actually called with;
#   * every numerical gradient agrees with a hand-derived exact one across six
#     surfaces and five points, and the worst error is reported, not hidden;
#   * of 360 bearings measured directly, the one that climbs fastest is the
#     gradient's -- to within the half-degree the sampling grid allows, and
#     the winning rate is |grad| times the cosine of that gap, to nine places;
#   * the gradient is perpendicular to an exactly parametrised contour, with
#     the dot product shrinking tenfold for each tenfold smaller step, which
#     is what "it goes to zero" looks like when every step is finite;
#   * a plane's gradient is the same vector everywhere -- until the function
#     value gets large, where roundoff eats it exactly as eps*|f|/2h predicts;
#   * three surfaces have the identical zero gradient at the origin and are a
#     minimum, a maximum and a saddle;
#   * the central difference's error on a cubic is exactly h squared, and the
#     total error is U-shaped in h with its trough at 1e-5;
#   * numpy.gradient differences a sampled array and is first-order at the
#     boundary by default, which is a different job from ours;
#   * nothing is left behind on disk.
#
# Everything after the one-time install runs offline. Nothing binds a port,
# nothing writes outside the lab, nothing needs a key. Deterministic,
# non-interactive, exits 0 only if every check passes.
set -u

export PYTHONDONTWRITEBYTECODE=1

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

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

failures=0
checks=0

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

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

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

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

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

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

echo "Day 109 — Which Way Is Uphill?"
echo

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

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

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

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

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

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

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

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

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

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

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

import surfaces as S
from gradients import (
    angle_degrees,
    angular_gap_degrees,
    contour_chord,
    directional_derivative,
    directional_derivative_direct,
    forward_partial,
    gradient,
    magnitude,
    partial,
    sweep_directions,
    unit,
)

def show(vector, places=6):
    return "[" + ",".join(f"{v:.{places}f}" for v in vector) + "]"

# -- partial derivatives ---------------------------------------------------
seen = []
def spy(p):
    seen.append(tuple(float(v) for v in p))
    return S.product(p)
partial(spy, (2.0, 5.0), 0)
print("evaluations_per_partial", len(seen))
print("frozen_coordinate_values", "|".join(str(v) for v in sorted({p[1] for p in seen})))
print("moved_coordinate_values", "|".join(f"{v:.5f}" for v in sorted(p[0] for p in seen)))

untouched = np.array([1.0, 2.0])
partial(S.bowl, untouched, 0)
print("point_unmutated", untouched.tolist() == [1.0, 2.0])
print("partial_is_a_plain_float", type(partial(S.bowl, (1.0, 1.0), 0)).__name__)

print("bowl_dfdx_at_2_1", round(partial(S.bowl, (2.0, 1.0), 0), 8))
print("bowl_dfdy_at_2_1", round(partial(S.bowl, (2.0, 1.0), 1), 8))
print("product_dfdx_at_1_0", round(partial(S.product, (1.0, 0.0), 0), 8))
print("product_dfdy_at_1_0", round(partial(S.product, (1.0, 0.0), 1), 8))

# -- gradients against exact -----------------------------------------------
worst = 0.0
count = 0
for name, (f, exact_gradient, _e, _g) in S.SURFACES.items():
    for p in S.PROBE_POINTS:
        worst = max(worst, float(np.max(np.abs(gradient(f, p) - exact_gradient(p)))))
        count += 1
print("gradients_checked", count)
print("worst_gradient_error_under_tolerance", worst < S.GRADIENT_TOL)
print("worst_gradient_error", f"{worst:.3e}")
print("tolerance_headroom_over_ten", S.GRADIENT_TOL / worst > 10.0)
print("bowl_gradient_at_1_1", show(gradient(S.bowl, (1.0, 1.0))))
print("cubic_gradient_at_2_1", show(gradient(S.cubic, (2.0, 1.0))))
print("gradient_size_two_inputs", gradient(S.bowl, (1.0, 1.0)).size)
print("gradient_size_three_inputs", gradient(S.model_loss, S.START_PARAMS).size)
print("bowl_gradient_magnitude", f"{magnitude(gradient(S.bowl, (1.0, 1.0))):.9f}")

# -- directional derivatives -----------------------------------------------
print("rate_due_east", round(directional_derivative(S.bowl, (1.0, 1.0), (1.0, 0.0)), 8))
print("rate_due_north", round(directional_derivative(S.bowl, (1.0, 1.0), (0.0, 1.0)), 8))
print("rate_along_3_minus_1", round(directional_derivative(S.bowl, (1.0, 1.0), (3.0, -1.0)), 8))
short = directional_derivative(S.bowl, (1.0, 1.0), (1.0, 2.0))
long = directional_derivative(S.bowl, (1.0, 1.0), (1000.0, 2000.0))
print("arrow_length_irrelevant", abs(short - long) < 1e-9)
worst_route = 0.0
for d in ((1.0, 1.0), (-1.0, 2.0), (3.0, -1.0), (-2.0, -5.0), (7.0, 0.5)):
    worst_route = max(worst_route, abs(
        directional_derivative(S.bowl, (1.0, 1.0), d)
        - directional_derivative_direct(S.bowl, (1.0, 1.0), d)))
print("dot_route_matches_direct_route", worst_route < S.GRADIENT_TOL)

# -- steepest ascent -------------------------------------------------------
angles, rates = sweep_directions(S.bowl, (1.0, 1.0))
best = int(np.argmax(rates))
worst_i = int(np.argmin(rates))
bearing = angle_degrees(S.bowl_gradient((1.0, 1.0)))
gap = angular_gap_degrees(float(np.degrees(angles[best])), bearing)
print("sweep_size", len(rates))
print("sweep_best_bearing", f"{np.degrees(angles[best]):.1f}")
print("gradient_bearing", f"{bearing:.4f}")
print("sweep_gap_degrees", f"{gap:.4f}")
print("sweep_gap_within_tolerance", gap <= S.ANGLE_TOL_DEGREES)
steepness = magnitude(S.bowl_gradient((1.0, 1.0)))
print("no_direction_beats_the_gradient", float(np.max(rates)) <= steepness + S.GRADIENT_TOL)
ratio = float(np.max(rates)) / steepness
print("best_rate_over_magnitude", f"{ratio:.12f}")
print("cosine_of_the_gap", f"{float(np.cos(np.radians(gap))):.12f}")
print("ratio_equals_cosine", abs(ratio - float(np.cos(np.radians(gap)))) < 1e-9)
print("up_and_down_are_opposite", f"{angular_gap_degrees(float(np.degrees(angles[best])), float(np.degrees(angles[worst_i]))):.4f}")
print("down_rate_is_minus_up_rate", abs(float(np.max(rates)) + float(np.min(rates))) < 1e-6)

# -- perpendicular to the contour ------------------------------------------
drift = 0.0
for name, (f, contour, level, _t0) in S.CONTOURS.items():
    for t in np.linspace(0.3, 2.6, 8):
        drift = max(drift, abs(f(contour(level, t)) - level))
print("contours_hold_f_constant", drift < 1e-12)
print("contour_drift", f"{drift:.3e}")

worst_dot = 0.0
for name, (f, contour, level, _t0) in S.CONTOURS.items():
    for t in (0.4, 0.9, 1.4, 1.9):
        chord, p, _q, _fp, _fq = contour_chord(f, contour, level, t, S.CONTOUR_DELTA)
        worst_dot = max(worst_dot, abs(float(np.dot(unit(gradient(f, p)), chord))))
print("perpendicular_within_tolerance", worst_dot < S.CONTOUR_DOT_TOL)
print("worst_contour_dot", f"{worst_dot:.3e}")

f, contour, level, t0 = S.CONTOURS["bowl"]
ratios = []
previous = None
for k in (2, 3, 4, 5, 6):
    chord, p, _q, _fp, _fq = contour_chord(f, contour, level, t0, 10.0 ** -k)
    dot = abs(float(np.dot(unit(gradient(f, p)), chord)))
    if previous is not None:
        ratios.append(previous / dot)
    previous = dot
print("dot_product_is_first_order_in_delta", all(9.0 < r < 11.0 for r in ratios))
print("dot_shrink_ratios", "|".join(f"{r:.3f}" for r in ratios))

a = np.sqrt(4.0)
b = np.sqrt(4.0 / 3.0)
tangent_worst = 0.0
for t in (0.0, 0.4, 0.9, 1.4, 1.9, 2.7):
    p = S.bowl_contour(4.0, t)
    tangent = np.array([-a * np.sin(t), b * np.cos(t)])
    tangent_worst = max(tangent_worst, abs(float(np.dot(tangent, S.bowl_gradient(p)))))
print("exact_tangent_dots_to_zero", tangent_worst < 1e-14)

# -- plane and bowl --------------------------------------------------------
seen_planes = [gradient(S.plane, p) for p in ((0.0, 0.0), (1.0, 1.0), (-40.0, 17.5))]
print("plane_gradient", show(seen_planes[0]))
print("plane_gradient_never_varies",
      float(np.max(np.abs(np.array(seen_planes) - seen_planes[0]))) < S.GRADIENT_TOL)
eps = float(np.finfo(float).eps)
far = (1000.0, -1000.0)
far_error = abs(partial(S.plane, far, 0) - 3.0)
predicted = eps * abs(S.plane(far)) / (2.0 * S.H_DEFAULT)
print("far_from_home_error_exceeds_the_labs_tolerance", far_error > S.GRADIENT_TOL)
print("far_from_home_error_matches_the_roundoff_bound",
      predicted / 100.0 < far_error < 3.0 * predicted)
print("far_error", f"{far_error:.3e}")
print("predicted_roundoff", f"{predicted:.3e}")

outward = all(float(np.dot(unit(gradient(S.bowl, p)), unit(np.array(p)))) > 0.0
              for p in ((0.5, 0.5), (1.0, 1.0), (2.0, 2.0), (-3.0, 1.0)))
print("bowl_gradient_points_away_from_the_minimum", outward)
lengths = [magnitude(gradient(S.bowl, (r, r))) for r in (0.5, 1.0, 2.0, 4.0)]
print("bowl_gradient_grows_with_distance", lengths == sorted(lengths))

# -- zero gradients --------------------------------------------------------
for name, kind, _why in S.STATIONARY_AT_ORIGIN:
    g = gradient(S.SURFACES[name][0], (0.0, 0.0))
    print(f"zero_gradient_{name}", f"{magnitude(g):.3e}")
print("saddle_walking_east", S.saddle((0.5, 0.0)))
print("saddle_walking_north", S.saddle((0.0, 0.5)))
print("saddle_on_the_diagonal", S.saddle((0.5, 0.5)))

# -- step size -------------------------------------------------------------
exact_dx = float(S.cubic_gradient((2.0, 1.0))[0])
print("cubic_exact_dfdx", exact_dx)
for k in (1, 2, 3):
    h = 10.0 ** -k
    err = partial(S.cubic, (2.0, 1.0), 0, h) - exact_dx
    print(f"cubic_error_is_h_squared_at_1e-{k}", abs(err - h * h) / (h * h) < 1e-5)
central = {k: abs(partial(S.cubic, (2.0, 1.0), 0, 10.0 ** -k) - exact_dx) for k in range(15)}
forward = {k: abs(forward_partial(S.cubic, (2.0, 1.0), 0, 10.0 ** -k) - exact_dx) for k in range(15)}
print("best_h_central", f"1e-{min(central, key=central.get):02d}")
print("best_h_forward", f"1e-{min(forward, key=forward.get):02d}")
print("central_beats_forward_at_default", forward[5] > 1000 * central[5])
print("tiny_h_is_worse_than_moderate_h", central[14] > central[1])
print("central_error_at_default", f"{central[5]:.3e}")
print("forward_error_at_default", f"{forward[5]:.3e}")

# -- numpy.gradient --------------------------------------------------------
xs = np.linspace(0.0, 4.0, 9)
ys = np.linspace(0.0, 4.0, 9)
X, Y = np.meshgrid(xs, ys, indexing="ij")
spacing = float(xs[1] - xs[0])
gx, gy = np.gradient(X * X + 3.0 * Y * Y, xs, ys)
print("npgradient_interior", f"[{gx[2, 2]:.6f},{gy[2, 2]:.6f}]")
print("npgradient_corner_default", f"[{gx[0, 0]:.6f},{gy[0, 0]:.6f}]")
gx2, gy2 = np.gradient(X * X + 3.0 * Y * Y, xs, ys, edge_order=2)
print("npgradient_corner_edge_order_2", f"[{gx2[0, 0]:.6f},{gy2[0, 0]:.6f}]")
cgx, _cgy = np.gradient(X ** 3 + X * Y * Y, xs, ys, edge_order=2)
cubic_exact = 3.0 * xs[4] ** 2 + ys[4] ** 2
print("npgradient_cubic_error_is_spacing_squared",
      abs(abs(cgx[4, 4] - cubic_exact) - spacing ** 2) < 1e-12)
print("npgradient_returns_a_field", gx.shape == (9, 9))
print("our_gradient_returns_a_vector", gradient(S.bowl, (1.0, 1.0)).shape == (2,))

# -- the model -------------------------------------------------------------
print("model_loss", S.model_loss(S.START_PARAMS))
print("model_gradient", show(S.model_loss_gradient(S.START_PARAMS), 1))
calls = []
def counted(p):
    calls.append(1)
    return S.model_loss(p)
gradient(counted, S.START_PARAMS)
print("evaluations_for_a_three_parameter_gradient", len(calls))
before = S.model_loss(S.START_PARAMS)
g = S.model_loss_gradient(S.START_PARAMS)
print("small_step_against_the_gradient_helps",
      S.model_loss(np.array(S.START_PARAMS) - 0.01 * g) < before)
print("small_step_along_the_gradient_hurts",
      S.model_loss(np.array(S.START_PARAMS) + 0.01 * g) > before)
print("too_large_a_step_overshoots",
      S.model_loss(np.array(S.START_PARAMS) - 0.2 * g) > before)
PY
)"

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

check_eq "a central difference evaluates f exactly twice" "2" "$(get evaluations_per_partial)"
check_eq "and the coordinate being held fixed never moves" "5.0" "$(get frozen_coordinate_values)"
check_eq "while the chosen coordinate moves h each way" \
  "1.99999|2.00001" "$(get moved_coordinate_values)"
check_eq "the caller's point is not mutated" "True" "$(get point_unmutated)"
check_eq "a partial returns a plain float, not a numpy scalar" \
  "float" "$(get partial_is_a_plain_float)"

check_eq "df/dx of x^2 + 3y^2 at (2, 1) is 4" "4.0" "$(get bowl_dfdx_at_2_1)"
check_eq "df/dy of the same is 6" "6.0" "$(get bowl_dfdy_at_2_1)"
check_eq "df/dx of xy at (1, 0) is 0" "0.0" "$(get product_dfdx_at_1_0)"
check_eq "df/dy of xy at the SAME point is 1, so the surface is not flat there" \
  "1.0" "$(get product_dfdy_at_1_0)"

check_eq "thirty gradients were checked against hand-derived exact ones" \
  "30" "$(get gradients_checked)"
check_eq "and every one is inside the stated tolerance" \
  "True" "$(get worst_gradient_error_under_tolerance)"
check_eq "with at least tenfold headroom rather than scraping past" \
  "True" "$(get tolerance_headroom_over_ten)"
echo "  (worst single gradient error on this run: $(get worst_gradient_error) -- reported, not asserted)"
check_eq "the bowl's gradient at (1, 1) is (2, 6)" \
  "[2.000000,6.000000]" "$(get bowl_gradient_at_1_1)"
check_eq "the cubic's gradient at (2, 1) is (13, 4)" \
  "[13.000000,4.000000]" "$(get cubic_gradient_at_2_1)"
check_eq "a two-input function has a two-component gradient" \
  "2" "$(get gradient_size_two_inputs)"
check_eq "and a three-input one has three" "3" "$(get gradient_size_three_inputs)"
check_eq "the gradient's length at (1, 1) is sqrt(40)" \
  "6.324555320" "$(get bowl_gradient_magnitude)"

check_eq "walking due east gives back the x partial" "2.0" "$(get rate_due_east)"
check_eq "walking due north gives back the y partial" "6.0" "$(get rate_due_north)"
check_eq "walking along (3, -1) gives exactly zero" "0.0" "$(get rate_along_3_minus_1)"
check_eq "a longer direction arrow does not give a bigger answer" \
  "True" "$(get arrow_length_irrelevant)"
check_eq "dotting with the gradient agrees with measuring along the direction" \
  "True" "$(get dot_route_matches_direct_route)"

check_eq "360 bearings were measured directly" "360" "$(get sweep_size)"
check_eq "and the fastest climb is at bearing 72" "72.0" "$(get sweep_best_bearing)"
check_eq "which is the gradient's own bearing to within a sampling step" \
  "71.5651" "$(get gradient_bearing)"
# Section 6 re-runs this script with D109_SELF_TEST=1, which swaps ONE
# expectation below for a deliberately wrong one. That is how the harness
# proves it can fail rather than merely asserting that it could.
expected_bearing="71.5651"
if [ -n "${D109_SELF_TEST:-}" ]; then
  expected_bearing="45.0000"   # the naive belief that on a bowl the gradient
                               # points straight away from the minimum
fi
check_eq "the gradient bearing on the bowl at (1, 1) is not the 45 degrees of the straight-back direction" \
  "${expected_bearing}" "$(get gradient_bearing)"
check_eq "the sampling gap is under the stated one-degree tolerance" \
  "True" "$(get sweep_gap_within_tolerance)"
check_eq "no direction anywhere beats the gradient's own magnitude" \
  "True" "$(get no_direction_beats_the_gradient)"
check_eq "the winning rate over the magnitude equals the cosine of the gap" \
  "True" "$(get ratio_equals_cosine)"
check_eq "the steepest descent is exactly 180 degrees round" \
  "180.0000" "$(get up_and_down_are_opposite)"
check_eq "and its rate is the negative of the steepest ascent" \
  "True" "$(get down_rate_is_minus_up_rate)"

check_eq "each parametrised contour really does hold f constant" \
  "True" "$(get contours_hold_f_constant)"
check_eq "the gradient is perpendicular to every contour tested" \
  "True" "$(get perpendicular_within_tolerance)"
echo "  (worst contour dot product on this run: $(get worst_contour_dot), tolerance 1e-04)"
check_eq "and the dot product shrinks tenfold for a tenfold smaller step" \
  "True" "$(get dot_product_is_first_order_in_delta)"
check_eq "the EXACT tangent and the EXACT gradient dot to zero, no tolerance needed" \
  "True" "$(get exact_tangent_dots_to_zero)"

check_eq "a plane's gradient is (3, -2)" "[3.000000,-2.000000]" "$(get plane_gradient)"
check_eq "and it is the same vector at every point tested" \
  "True" "$(get plane_gradient_never_varies)"
check_eq "far from the origin the estimate breaks the lab's own tolerance" \
  "True" "$(get far_from_home_error_exceeds_the_labs_tolerance)"
check_eq "by an amount the roundoff bound eps|f|/2h predicts" \
  "True" "$(get far_from_home_error_matches_the_roundoff_bound)"
echo "  (measured $(get far_error) against a predicted $(get predicted_roundoff))"

check_eq "the bowl's gradient points away from its minimum, everywhere tested" \
  "True" "$(get bowl_gradient_points_away_from_the_minimum)"
check_eq "and gets longer the further out you stand" \
  "True" "$(get bowl_gradient_grows_with_distance)"

check_eq "the bowl has a zero gradient at the origin" "0.000e+00" "$(get zero_gradient_bowl)"
check_eq "so does the dome" "0.000e+00" "$(get zero_gradient_dome)"
check_eq "so does the saddle" "0.000e+00" "$(get zero_gradient_saddle)"
check_eq "yet walking east from the saddle goes UP" "0.25" "$(get saddle_walking_east)"
check_eq "and walking north from it goes DOWN" "-0.25" "$(get saddle_walking_north)"
check_eq "and along the diagonal nothing happens at all" \
  "0.0" "$(get saddle_on_the_diagonal)"

check_eq "the cubic's exact df/dx at (2, 1) is 13" "13.0" "$(get cubic_exact_dfdx)"
check_eq "the central difference overshoots by exactly h^2 at h = 1e-1" \
  "True" "$(get cubic_error_is_h_squared_at_1e-1)"
check_eq "and at h = 1e-2" "True" "$(get cubic_error_is_h_squared_at_1e-2)"
check_eq "and at h = 1e-3" "True" "$(get cubic_error_is_h_squared_at_1e-3)"
check_eq "the best central step over 15 decades is 1e-05" "1e-05" "$(get best_h_central)"
check_eq "the best forward step is 1e-08, three decades away" "1e-08" "$(get best_h_forward)"
check_eq "at the default step central beats forward by over a thousandfold" \
  "True" "$(get central_beats_forward_at_default)"
check_eq "and a step of 1e-14 is worse than one of 1e-01" \
  "True" "$(get tiny_h_is_worse_than_moderate_h)"
echo "  (central $(get central_error_at_default) against forward $(get forward_error_at_default) at h = 1e-5)"

check_eq "numpy.gradient is exact in the interior of a sampled quadratic" \
  "[2.000000,6.000000]" "$(get npgradient_interior)"
check_eq "but first-order at the corner by default, giving (0.5, 1.5) where (0, 0) is right" \
  "[0.500000,1.500000]" "$(get npgradient_corner_default)"
check_eq "edge_order=2 fixes that corner exactly" \
  "[0.000000,0.000000]" "$(get npgradient_corner_edge_order_2)"
check_eq "on a cubic its error is the GRID spacing squared, which you cannot choose" \
  "True" "$(get npgradient_cubic_error_is_spacing_squared)"
check_eq "numpy.gradient returns a field over the whole array" \
  "True" "$(get npgradient_returns_a_field)"
check_eq "ours returns one vector at one point" \
  "True" "$(get our_gradient_returns_a_vector)"

check_eq "the three-parameter loss is 22.5" "22.5" "$(get model_loss)"
check_eq "and its gradient is the three whole numbers (-17, -18, -8)" \
  "[-17.0,-18.0,-8.0]" "$(get model_gradient)"
check_eq "which cost six evaluations of the loss: two per parameter" \
  "6" "$(get evaluations_for_a_three_parameter_gradient)"
check_eq "a small step AGAINST the gradient reduces the loss" \
  "True" "$(get small_step_against_the_gradient_helps)"
check_eq "a small step ALONG it increases the loss" \
  "True" "$(get small_step_along_the_gradient_hurts)"
check_eq "and too large a step overshoots to worse than the start" \
  "True" "$(get too_large_a_step_overshoots)"

# --------------------------------------------------------------------------
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 a wrong one -- 45 degrees, which is where the gradient WOULD point if a
# bowl's uphill direction were simply "straight away from the bottom" -- and
# asserts that the re-run reports the failure and exits non-zero. If this
# section passes, section 5 is not decorative.
if [ -z "${D109_SELF_TEST:-}" ]; then
  self_out="$(D109_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: the gradient bearing on the bowl at (1, 1) is not the 45 degrees"*)
      check "the failing check is named in the output with both values" "yes" ;;
    *) check "the failing check is named in the output with both values" "no" ;;
  esac
  case "${self_out}" in
    *", 1 failure(s)."*)
      check "the summary line counts exactly one failure" "yes" ;;
    *) check "the summary line counts exactly one failure" "no" ;;
  esac
else
  echo "  (self-test run: section 6 does not recurse)"
fi

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

# `.venv` is pruned from both searches below. The virtual environment ships
# NumPy's and pytest's own precompiled bytecode -- hundreds of __pycache__
# directories that came with the packages and have nothing to do with whether
# THIS lab tidied up after itself. Searching them would report a failure the
# reader cannot fix and did not cause. Everything the lab itself writes lives
# outside `.venv`, which is exactly what these two checks look at.

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 by the lab's own code" "no"
else
  check "no __pycache__ directory left by the lab's own code" "yes"
fi

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

# And a check on the check. If a future edit dropped the `-name '.venv'
# -prune` from either search above, the two checks would start reporting a
# failure caused entirely by NumPy's own shipped bytecode. This proves the
# prune is doing its job: when a `.venv` exists and contains __pycache__
# directories of its own, the pruned search must still find nothing.
if [ -d "${lab_dir}/.venv" ]; then
  inside="$(find "${lab_dir}/.venv" -type d -name '__pycache__' -print -quit 2>/dev/null)"
  if [ -n "${inside}" ]; then
    outside="$(find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -print -quit 2>/dev/null)"
    if [ -z "${outside}" ]; then
      check "the .venv prune works: the environment's own bytecode is not counted against the lab" "yes"
    else
      check "the .venv prune works: the environment's own bytecode is not counted against the lab" "no"
    fi
  else
    check "the lab-local .venv exists and holds no bytecode caches of its own" "yes"
  fi
else
  check "no lab-local .venv on this run, so there is nothing to prune" "yes"
fi

# `.venv` must never be reported as a stray file either. It is created by the
# documented setup commands in the README; a suite that then complained about
# its existence would be telling the reader off for following instructions.
if [ -d "${lab_dir}/.venv" ]; then
  check ".venv is treated as expected, never as something left behind" "yes"
else
  check ".venv is treated as expected, never as something left behind" "yes"
fi

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

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