Hands-on lab — Day 108: Derivatives: Rates of Change

"""A rate is a difference divided by the interval it happened over.

Run from inside `examples/`:

    ../.venv/bin/python3 01_average_rate_of_change.py

Nothing in this script is calculus yet. It is a stopwatch and a tape measure.
"""

from __future__ import annotations

import dataset as D
from derivatives import average_rate

print("Day 108 / 01 — average rate of change")
print()

# --------------------------------------------------------------------------
print("1. A car, timed once a second")
# --------------------------------------------------------------------------
print()
print("   The distances are invented. They are 4 * t**2 metres, so the car is")
print("   speeding up steadily and every number below can be checked by hand.")
print()
print("     t (s)   distance (m)")
for t, d in zip(D.CAR_TIMES_S, D.CAR_DISTANCE_M):
    print(f"     {t:5.1f}   {d:12.1f}")
print()

# --------------------------------------------------------------------------
print("2. Rise over run, over the whole trip")
# --------------------------------------------------------------------------
print()


def distance(t: float) -> float:
    """Distance in metres at time t seconds. 4 * t**2, matching the table."""
    return 4.0 * t * t


whole = average_rate(distance, 0.0, 6.0)
print("     rise = 144.0 - 0.0 =", D.CAR_DISTANCE_M[-1] - D.CAR_DISTANCE_M[0], "metres")
print("     run  =   6.0 - 0.0 =", D.CAR_TIMES_S[-1] - D.CAR_TIMES_S[0], "seconds")
print(f"     average speed = rise / run = {whole} m/s")
assert whole == D.CAR_AVERAGE_SPEED_WHOLE_TRIP
print()
print("   24 metres per second, averaged over six seconds. Note what that")
print("   number does NOT say: the car was never travelling at 24 m/s for the")
print("   whole trip. It started at rest and finished much faster.")
print()

# --------------------------------------------------------------------------
print("3. The same question, asked of one second at a time")
# --------------------------------------------------------------------------
print()
print("     interval      rise (m)   run (s)   average speed (m/s)")
per_second = []
for i in range(len(D.CAR_TIMES_S) - 1):
    a, b = D.CAR_TIMES_S[i], D.CAR_TIMES_S[i + 1]
    rise = D.CAR_DISTANCE_M[i + 1] - D.CAR_DISTANCE_M[i]
    speed = average_rate(distance, a, b)
    per_second.append(speed)
    print(f"     [{a:.0f}, {b:.0f}]        {rise:8.1f}   {b - a:7.1f}   {speed:19.1f}")
print()
assert per_second == [4.0, 12.0, 20.0, 28.0, 36.0, 44.0]
assert per_second[3] == D.CAR_AVERAGE_SPEED_SECOND_FOUR
print("   Six different answers to 'how fast was the car'. All six are correct.")
print("   They are answers to six different questions.")
print()

# --------------------------------------------------------------------------
print("4. The question a speedometer answers")
# --------------------------------------------------------------------------
print()
print("   A speedometer does not show an average over an interval. It shows a")
print("   number at an INSTANT. Ask for the average speed over an interval of")
print("   no width at all and the arithmetic refuses:")
print()
try:
    average_rate(distance, 3.0, 3.0)
except ZeroDivisionError as exc:
    print("     ZeroDivisionError:", str(exc).split(";")[0])
else:  # pragma: no cover - the call above always raises
    raise AssertionError("average_rate over a zero-width interval must refuse")
print()
print("   Rise zero, run zero, and 0/0 is not a number. That refusal is the")
print("   whole problem, and the derivative is the machine that gets round it:")
print("   instead of asking for the rate over no interval, ask for the rate")
print("   over intervals that get smaller and smaller, and see whether the")
print("   answers settle on something. Script 02 does exactly that.")
print()
print("   For this car the settled answer at t = 3 is", D.CAR_INSTANT_SPEED_AT_3, "m/s,")
print("   which sits between the 20 m/s averaged over second three and the")
print("   28 m/s averaged over second four -- as it must.")
assert per_second[2] < D.CAR_INSTANT_SPEED_AT_3 < per_second[3]
print()

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

"""Shrink the interval and watch the answer settle. That settling is the limit.

Run from inside `examples/`:

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

from __future__ import annotations

import dataset as D
from derivatives import average_rate, shrinking_slopes, tangent_at

print("Day 108 / 02 — the shrinking interval")
print()

# --------------------------------------------------------------------------
print("1. The car again, at t = 3, over intervals that get smaller")
# --------------------------------------------------------------------------
print()


def distance(t: float) -> float:
    return 4.0 * t * t


print("     width h    interval          average speed (m/s)")
for h in [1.0, 0.5, 0.1, 0.01, 0.001, 0.0001]:
    speed = average_rate(distance, 3.0, 3.0 + h)
    print(f"     {h:<10.4f} [3, {3.0 + h:<8.4f}]   {speed:19.6f}")
print()
print("   The numbers are 24 + 4h, and you can see them heading for 24 without")
print("   being told. Nobody computed a limit here; the sequence was watched.")
print()

# --------------------------------------------------------------------------
print("2. The same thing on f(x) = x**2 at x = 3, with the algebra shown")
# --------------------------------------------------------------------------
print()
print("     f(3 + h) - f(3)     (3 + h)**2 - 9      9 + 6h + h**2 - 9")
print("     ----------------  =  --------------  =  ------------------  =  6 + h")
print("            h                   h                    h")
print()
slopes = shrinking_slopes(D.square, D.SETTLE_POINT, D.SETTLE_WIDTHS)
print("     width h        secant slope     exactly 6 + h?   distance from 6")
for h, slope, expected in zip(D.SETTLE_WIDTHS, slopes, D.SETTLE_EXPECTED_SLOPES):
    matches = abs(slope - expected) < D.EXACT_TOL
    print(f"     {h:<12.4f}   {slope:<16.12f} {str(matches):<16} {abs(slope - 6.0):.4f}")
    assert matches, (h, slope, expected)
print()
print("   Two facts sit in that table and they are not the same fact.")
print()
print("   The first is that the slope over [3, 3 + h] is 6 + h EXACTLY, for")
print("   every h, with no approximation anywhere. That is algebra.")
print()
print("   The second is that as h shrinks, 6 + h gets arbitrarily close to 6.")
print("   That is the limit, and it is why the derivative of x**2 at 3 is 6.")
print()
print("   Notice what the algebra did that the arithmetic could not: it")
print("   cancelled the h in the denominator BEFORE h was allowed to reach")
print("   zero. At h = 0 the fraction is 0/0 and means nothing; the simplified")
print("   form 6 + h means something at every h including zero.")
print()

# --------------------------------------------------------------------------
print("3. The floating-point version of the same sequence is not exact")
# --------------------------------------------------------------------------
print()
print("     h            computed slope          6 + h            gap")
for h, slope, expected in zip(D.SETTLE_WIDTHS, slopes, D.SETTLE_EXPECTED_SLOPES):
    print(f"     {h:<12.4f} {slope!r:<24} {expected:<16.4f} {abs(slope - expected):.3e}")
print()
print("   The gaps are around 1e-13, and they are not mistakes in the formula.")
print("   They are the arithmetic: 3.001**2 cannot be stored exactly in binary,")
print("   the subtraction loses some of the digits the two numbers had in")
print("   common, and dividing by 0.001 multiplies what is left by a thousand.")
print("   Script 05 measures that effect properly, because it is the reason a")
print("   smaller h eventually makes a numerical derivative WORSE.")
print()

# --------------------------------------------------------------------------
print("4. Secants approaching a tangent")
# --------------------------------------------------------------------------
print()
print("   Each row is a straight line through (3, 9) and one other point on")
print("   the curve. As the other point slides in, the line pivots.")
print()
def line_text(slope: float, intercept: float) -> str:
    """y = mx + c, written with the sign of c folded into the operator."""
    sign = "-" if intercept < 0 else "+"
    return f"y = {slope:.2f}x {sign} {abs(intercept):.2f}"


print("     h        second point           slope    line through (3, 9)")
for h in [2.0, 1.0, 0.5, 0.1]:
    x2 = 3.0 + h
    slope = average_rate(D.square, 3.0, x2)
    intercept = 9.0 - slope * 3.0
    print(f"     {h:<8.2f} ({x2:.2f}, {D.square(x2):7.4f})       {slope:6.2f}   {line_text(slope, intercept)}")
print()
tangent_slope, tangent_intercept = tangent_at(D.square, 3.0, 1e-5)
print(f"     tangent (h -> 0)                       {tangent_slope:6.2f}   "
      f"{line_text(round(tangent_slope, 2), round(tangent_intercept, 2))}")
assert abs(tangent_slope - 6.0) < D.CENTRAL_TOL
assert abs(tangent_intercept - (-9.0)) < 1e-8
print()
print("   The tangent is not 'the line that touches the curve once' -- plenty")
print("   of lines do that and are not tangents. It is the line the secants")
print("   approach, and its slope is the derivative. y = 6x - 9 touches the")
print("   parabola at x = 3 and matches its direction there.")
print()

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

"""The handful of rules you actually need, each one checked against a number.

Run from inside `examples/`:

    ../.venv/bin/python3 03_rules_checked_numerically.py

A rule you have only been told is a rule you cannot check. Every one below is
stated, then computed numerically at a point, then compared.
"""

from __future__ import annotations

import math

import dataset as D
from derivatives import central_difference

print("Day 108 / 03 — the rules, and the numbers that agree with them")
print()

# --------------------------------------------------------------------------
print("1. Notation, defined once")
# --------------------------------------------------------------------------
print()
print("   Three ways of writing the same object, all in current use:")
print()
print("     f'(x)     'f prime of x'. Lagrange's notation. Compact, and the")
print("               one to reach for when the input variable is obvious.")
print("     dy/dx     Leibniz's notation, read 'dee y by dee x'. It names the")
print("               two variables, which matters the moment there is more")
print("               than one input -- Day 109's whole subject.")
print("     Df(x)     Euler's operator notation. You will meet it in papers.")
print()
print("   dy/dx is not a fraction, although it is descended from one and")
print("   behaves like one often enough to be dangerous. Read it as: the rate")
print("   at which y changes with respect to x.")
print()

# --------------------------------------------------------------------------
print("2. The rules")
# --------------------------------------------------------------------------
print()
print("     constant           d/dx of c          =  0")
print("     power              d/dx of x**n       =  n * x**(n-1)")
print("     constant multiple  d/dx of c * f(x)   =  c * f'(x)")
print("     sum                d/dx of f(x)+g(x)  =  f'(x) + g'(x)")
print("     exponential        d/dx of e**x       =  e**x")
print("     logarithm          d/dx of ln(x)      =  1/x")
print()
print("   The first four are worth understanding. The constant rule says a")
print("   flat line has no slope. The power rule generalises the (3+h)**2")
print("   expansion from script 02: multiply out, cancel the h, and what")
print("   survives is n * x**(n-1). The constant-multiple rule says stretching")
print("   a graph vertically by 5 stretches every slope by 5. The sum rule says")
print("   rates add, which is why it is safe to differentiate a long expression")
print("   one term at a time.")
print()
print("   The last two are facts to know. They are not obvious and today does")
print("   not derive them.")
print()

# --------------------------------------------------------------------------
print("3. Every rule, checked at a point")
# --------------------------------------------------------------------------
print()
print(f"   Numerically, with the central difference at h = {D.COMPARE_WIDTH}.")
print("   'exact' is what the rule says. 'measured' is what the arithmetic")
print("   says without knowing the rule.")
print()
print("     rule                                        at x    exact          measured       error")
for (name, f, exact_derivative, x), expected in zip(D.RULE_CASES, D.RULE_EXPECTED):
    exact = exact_derivative(x)
    measured = central_difference(f, x, D.COMPARE_WIDTH)
    error = abs(measured - exact)
    print(f"     {name:<42}  {x:<6.2f}  {exact:<14.9f} {measured:<14.9f} {error:.2e}")
    assert exact == expected or abs(exact - expected) < 1e-15, (name, exact, expected)
    assert error < D.RULE_TOL, (name, error)
print()
print(f"   Every error is below {D.RULE_TOL:.0e}, and that ceiling was not chosen by")
print("   running the check and enlarging the number until it passed. The")
print("   central difference's truncation error is about h**2/6 times the third")
print("   derivative; the worst case here is x**5 at 1.5, where the third")
print("   derivative is 135, giving 1e-10/6 * 135 = 2.25e-9. The ceiling is")
print("   four times that. See dataset.py, which shows the arithmetic.")
print()

# --------------------------------------------------------------------------
print("4. Why e is the special one")
# --------------------------------------------------------------------------
print()
print("   Take b**x for various bases and measure its slope at x = 0.")
print("   The value at x = 0 is 1 for every base, so the slope is the only")
print("   thing that distinguishes them.")
print()
print("     base b        f(0)      f'(0) measured")
ratios = {}
for base in [2.0, 2.5, math.e, 3.0, 10.0]:
    def power_of(x: float, b: float = base) -> float:
        return b**x

    slope = central_difference(power_of, 0.0, D.COMPARE_WIDTH)
    ratios[base] = slope
    label = "e = 2.718..." if base == math.e else f"{base}"
    print(f"     {label:<13} {power_of(0.0):<9.1f} {slope:.9f}")
print()
assert ratios[2.0] < 1.0 < ratios[3.0]
assert abs(ratios[math.e] - 1.0) < D.CENTRAL_TOL
print("   For base 2 the slope at 0 is less than 1. For base 3 it is more than")
print("   1. Somewhere between them is a base whose slope at 0 is EXACTLY 1,")
print("   and that base is e. Measured here as", f"{ratios[math.e]:.12f}.")
print()
print("   That is what makes e**x its own derivative. Every b**x is")
print("   proportional to its own derivative -- the graph's steepness is always")
print("   proportional to its height -- and e is the base for which the")
print("   constant of proportionality is 1 rather than 0.693 or 1.099.")
print()
print("   Check the proportionality directly, at three different points:")
print()
print("     x       e**x            measured f'(x)   ratio f'(x) / f(x)")
for x in [0.0, 1.0, 2.5]:
    value = math.exp(x)
    slope = central_difference(math.exp, x, D.COMPARE_WIDTH)
    ratio = slope / value
    print(f"     {x:<7.1f} {value:<15.9f} {slope:<16.9f} {ratio:.12f}")
    assert abs(ratio - 1.0) < 1e-9, (x, ratio)
print()
print("   The ratio is 1 everywhere, which is the whole claim.")
print()

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

"""Two ways to estimate a slope from function values, and why one is far better.

Run from inside `examples/`:

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

from __future__ import annotations

import math

import dataset as D
from derivatives import (
    backward_difference,
    central_difference,
    forward_difference,
    numpy_gradient_slope,
    numpy_gradient_slope_from_coordinates,
)

print("Day 108 / 04 — forward, backward, central")
print()

# --------------------------------------------------------------------------
print("1. The three rules")
# --------------------------------------------------------------------------
print()
print("     forward    ( f(x + h) - f(x)     ) / h        one step ahead")
print("     backward   ( f(x)     - f(x - h) ) / h        one step behind")
print("     central    ( f(x + h) - f(x - h) ) / (2h)     straddling x")
print()
print("   The central rule is the average of the other two, and the averaging")
print("   is the entire trick. Forward leans one way off the tangent, backward")
print("   leans the other way by almost exactly the same amount, and adding")
print("   them cancels the leaning. It costs one function call more than")
print("   forward and none at all more than computing both one-sided rules.")
print()

# --------------------------------------------------------------------------
print("2. All three on f(x) = x**2 at x = 3, where the answer is 6")
# --------------------------------------------------------------------------
print()
print("     h         forward       backward      central       f-err     c-err")
for h in [1.0, 0.1, 0.01, 0.001]:
    fwd = forward_difference(D.square, 3.0, h)
    bwd = backward_difference(D.square, 3.0, h)
    cen = central_difference(D.square, 3.0, h)
    print(f"     {h:<9.4f} {fwd:<13.9f} {bwd:<13.9f} {cen:<13.9f} {abs(fwd - 6.0):<9.2e} {abs(cen - 6.0):.1e}")
    assert abs(fwd - (6.0 + h)) < 1e-9
    assert abs(bwd - (6.0 - h)) < 1e-9
    assert abs(cen - 6.0) < 1e-9
print()
print("   Forward is 6 + h. Backward is 6 - h. Their average is 6, and for a")
print("   parabola that is not an approximation -- it is exact, at every h.")
print("   The forward rule's error is proportional to h; the central rule's")
print("   error on a quadratic is zero because a quadratic has no third")
print("   derivative for it to trip over.")
print()

# --------------------------------------------------------------------------
print("3. On e**x at x = 1, where central is very good but not exact")
# --------------------------------------------------------------------------
print()
print(f"   The right answer is e = {math.e!r}")
print()
print("     h            forward error    central error    central is better by")
for h in [1e-1, 1e-2, 1e-3, 1e-4, 1e-5]:
    fwd_err = abs(forward_difference(D.exponential, 1.0, h) - math.e)
    cen_err = abs(central_difference(D.exponential, 1.0, h) - math.e)
    print(f"     {h:<12.0e} {fwd_err:<16.6e} {cen_err:<16.6e} {fwd_err / cen_err:>10.0f}x")
    assert cen_err < fwd_err
print()
print("   Read down the two error columns. Each time h drops by a factor of 10,")
print("   the forward error drops by about 10 and the central error by about")
print("   100. That is the difference between an error proportional to h and an")
print("   error proportional to h**2, and it compounds: at h = 1e-5 the central")
print("   rule is over two hundred thousand times more accurate for one extra")
print("   function call.")
print()
h = D.COMPARE_WIDTH
fwd_err = abs(forward_difference(D.exponential, 1.0, h) - math.e)
cen_err = abs(central_difference(D.exponential, 1.0, h) - math.e)
assert fwd_err < D.FORWARD_TOL
assert cen_err < D.CENTRAL_TOL
assert cen_err * 1000.0 < fwd_err
print(f"   At h = {h:.0e}: forward error {fwd_err:.6e}, central error {cen_err:.6e}.")
print(f"   Both are inside the tolerances derived in dataset.py "
      f"({D.FORWARD_TOL:.0e} and {D.CENTRAL_TOL:.0e}).")
print()

# --------------------------------------------------------------------------
print("4. Where the h**2 comes from, without any calculus")
# --------------------------------------------------------------------------
print()
print("   Taylor's expansion writes a smooth function near x as a polynomial:")
print()
print("     f(x + h) = f(x) + h*f'(x) + (h**2/2)*f''(x) + (h**3/6)*f'''(x) + ...")
print("     f(x - h) = f(x) - h*f'(x) + (h**2/2)*f''(x) - (h**3/6)*f'''(x) + ...")
print()
print("   Subtract the first from f(x) and divide by h and the f'' term")
print("   survives, multiplied by h/2. That is the forward rule's error.")
print()
print("   Subtract the SECOND from the first and the f'' terms cancel -- they")
print("   have the same sign in both lines. Divide by 2h and the first survivor")
print("   is the f''' term, multiplied by h**2/6. That is the central rule's")
print("   error, and it is why halving h quarters it.")
print()
print("   Check the prediction against the measurement on e**x at x = 1, where")
print("   every derivative equals e:")
print()
print("     h         predicted h**2/6 * e    measured central error   ratio")
for h_ in [1e-2, 1e-3, 1e-4]:
    predicted = (h_ * h_ / 6.0) * math.e
    measured = abs(central_difference(D.exponential, 1.0, h_) - math.e)
    print(f"     {h_:<9.0e} {predicted:<23.6e} {measured:<24.6e} {measured / predicted:.4f}")
    assert 0.99 < measured / predicted < 1.01, (h_, measured / predicted)
print()
print("   Within one percent, three times over. The formula is not folklore.")
print()

# --------------------------------------------------------------------------
print("5. The same job, handed to NumPy")
# --------------------------------------------------------------------------
print()
print("   numpy.gradient differentiates SAMPLES, not functions: you give it")
print("   values you already have and it returns a slope estimate at each one.")
print("   Interior points get the central rule; the two ends have nothing on")
print("   one side and fall back to a one-sided rule.")
print()
mine = central_difference(D.exponential, 1.0, D.COMPARE_WIDTH)
theirs = numpy_gradient_slope(D.exponential, 1.0, D.COMPARE_WIDTH)
print(f"     central_difference        {mine!r}")
print(f"     np.gradient (interior)    {theirs!r}")
print(f"     identical to the last bit {mine == theirs}")
assert mine == theirs
print()
coords = numpy_gradient_slope_from_coordinates(D.exponential, 1.0, D.COMPARE_WIDTH)
print("   One honest wrinkle, found while building this lab. Pass np.gradient a")
print("   scalar spacing and it is bit-for-bit our central difference. Pass it")
print("   an ARRAY of coordinates -- the same evenly spaced points -- and it")
print("   uses its general unevenly-spaced formula instead:")
print()
print(f"     np.gradient(ys, h)        {theirs!r}")
print(f"     np.gradient(ys, xs)       {coords!r}")
print(f"     they differ by            {abs(coords - theirs):.3e}")
assert coords != theirs
assert abs(coords - theirs) < 1e-10
print()
print("   Both are correct. Neither is the exact answer. 'The same formula' is")
print("   a claim about the mathematics, and the mathematics does not fix the")
print("   order the additions happen in.")
print()

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

"""Making h smaller makes the answer better -- until it makes it much worse.

Run from inside `examples/`:

    ../.venv/bin/python3 05_the_u_shaped_error.py

This is the script the lab exists for. It contradicts the obvious intuition,
and it contradicts it with a measurement rather than an argument.
"""

from __future__ import annotations

import math

import dataset as D
from derivatives import (
    best_step,
    central_difference,
    error_curve,
    forward_difference,
    is_u_shaped,
    numpy_error_curve,
)

print("Day 108 / 05 — the U-shaped error curve")
print()

# --------------------------------------------------------------------------
print("1. The intuition, and why it is wrong")
# --------------------------------------------------------------------------
print()
print("   The derivative is the limit of the difference quotient as h goes to")
print("   zero. So a smaller h should give a better answer, and h = 1e-300")
print("   should give a nearly perfect one.")
print()
print("   It gives 0.0. Here it is:")
print()
tiny = forward_difference(D.exponential, 1.0, 1e-300)
print(f"     forward_difference(exp, 1.0, 1e-300)  ->  {tiny!r}")
print(f"     the right answer is                       {math.e!r}")
assert tiny == 0.0
print()
print("   Not slightly wrong. Not wrong in the eighth decimal place. Zero,")
print("   with total confidence and no warning of any kind.")
print()
print("   The reason: exp(1 + 1e-300) and exp(1) are the same float64. There is")
print("   no float between them, so their difference is exactly 0.0, and 0.0")
print("   divided by anything is 0.0. The subtraction destroyed every digit the")
print("   two numbers had in common -- which was all of them.")
print()
same = math.exp(1.0 + 1e-300) == math.exp(1.0)
print(f"     exp(1 + 1e-300) == exp(1)  ->  {same}")
assert same
print()

# --------------------------------------------------------------------------
print("2. Two errors, pulling in opposite directions")
# --------------------------------------------------------------------------
print()
print("   TRUNCATION error comes from the mathematics. The formula is the")
print("   limit's approximation at a finite h, and it is wrong by roughly")
print("     forward:  (h/2)    * f''(x)")
print("     central:  (h**2/6) * f'''(x)")
print("   Both SHRINK as h shrinks. This is the term everybody knows about.")
print()
print("   ROUNDING error comes from the arithmetic. f(x+h) and f(x-h) are each")
print("   stored to about 1e-16 relative precision. Subtracting two nearly")
print("   equal numbers throws away their leading digits and leaves the noise;")
print("   dividing by a tiny h then multiplies that noise by 1/h:")
print("     either rule:  about  EPSILON * |f(x)| / h")
print("   This term GROWS as h shrinks. This is the term that surprises people.")
print()
print(f"     EPSILON (float64) = {D.EPSILON!r}")
print()
print("   Add a term that shrinks to a term that grows and you get a U. The")
print("   bottom of the U is the best h there is, and it is nowhere near zero.")
print()

# --------------------------------------------------------------------------
print("3. The measurement: 27 step sizes from 1e-1 down to 1e-14")
# --------------------------------------------------------------------------
print()
print("   f(x) = e**x at x = 1. The exact slope is e, so the error is knowable.")
print()
forward_errors = error_curve(D.exponential, D.U_POINT, D.U_EXACT_SLOPE, D.U_WIDTHS, forward_difference)
central_errors = error_curve(D.exponential, D.U_POINT, D.U_EXACT_SLOPE, D.U_WIDTHS, central_difference)

print("     h            forward error   central error   central error, log-log")
scale = 3.0
for h, fe, ce in zip(D.U_WIDTHS, forward_errors, central_errors):
    exponent = math.log10(ce) if ce > 0 else -18.0
    bar_length = max(0, int(round((exponent + 12.0) * scale)))
    print(f"     {h:<12.3e} {fe:<15.6e} {ce:<15.6e} {'#' * bar_length}")
print()
print("   Read the bars. They shorten as h shrinks -- the error is falling --")
print("   reach their shortest in the middle of the table, and then lengthen")
print("   again all the way to the bottom. That is the U, drawn sideways.")
print()

# --------------------------------------------------------------------------
print("4. Where the bottom is")
# --------------------------------------------------------------------------
print()
best_forward_h, best_forward_error = best_step(D.U_WIDTHS, forward_errors)
best_central_h, best_central_error = best_step(D.U_WIDTHS, central_errors)
print(f"     forward:  best h = {best_forward_h:.3e}   error there = {best_forward_error:.6e}")
print(f"     central:  best h = {best_central_h:.3e}   error there = {best_central_error:.6e}")
print()
assert is_u_shaped(forward_errors)
assert is_u_shaped(central_errors)
assert best_central_error < best_forward_error
assert forward_errors[0] > 10.0 * best_forward_error
assert forward_errors[-1] > 10.0 * best_forward_error
assert central_errors[0] > 10.0 * best_central_error
assert central_errors[-1] > 10.0 * best_central_error
print("   Both curves are U-shaped by the test in derivatives.py: the minimum")
print("   is in the interior, and both ends are more than ten times worse than")
print("   the middle. Both facts are asserted, not eyeballed.")
print()
print("   Balancing the two error terms predicts where the bottom should be.")
print("   Setting (h/2)*e equal to EPSILON*e/h and solving gives h about")
print("   sqrt(2*EPSILON) for the forward rule; setting (h**2/6)*e equal to")
print("   EPSILON*e/h gives h about (3*EPSILON)**(1/3) for the central rule.")
print()
predicted_forward = math.sqrt(2.0 * D.EPSILON)
predicted_central = (3.0 * D.EPSILON) ** (1.0 / 3.0)
print(f"     forward:  predicted {predicted_forward:.3e}   measured {best_forward_h:.3e}")
print(f"     central:  predicted {predicted_central:.3e}   measured {best_central_h:.3e}")
print()
assert 0.1 < best_forward_h / predicted_forward < 10.0
assert 0.1 < best_central_h / predicted_central < 10.0
print("   Both measurements land within a factor of ten of the prediction, and")
print("   a factor of ten is the right expectation: the grid here has three")
print("   steps per decade, the constants in the two error terms were dropped,")
print("   and the rounding term is a random walk rather than a smooth curve.")
print("   The lab asserts the order of magnitude and reports the rest.")
print()
print("   The practical rule that falls out, for float64:")
print()
print("     forward difference   h around 1e-8")
print("     central difference   h around 1e-5 to 1e-6")
print()
print("   And the practical warning: h = 1e-12 is not a careful choice. It is")
print(f"   a worse answer than h = 1e-3, by a factor of "
      f"{central_errors[D.U_WIDTHS.index(1e-12)] / central_errors[D.U_WIDTHS.index(1e-3)]:,.0f} "
      "on this run.")
print()

# --------------------------------------------------------------------------
print("5. The noise at the bottom is real noise")
# --------------------------------------------------------------------------
print()
print("   The error does not fall smoothly to a point and rise smoothly away.")
print("   Around the minimum it jitters, because the rounding term depends on")
print("   exactly which bits happen to survive the subtraction at that h:")
print()
print("     h            central error")
for h, ce in zip(D.U_WIDTHS, central_errors):
    if 1e-8 <= h <= 1e-5:
        print(f"     {h:<12.3e} {ce:.6e}")
print()
print("   That is why the lab's U-shape test asks about the ends against the")
print("   middle rather than for a monotone descent. A test demanding smooth")
print("   monotonicity here would be a test demanding something untrue.")
print()

# --------------------------------------------------------------------------
print("6. The array version, for the plot")
# --------------------------------------------------------------------------
print()
array_errors = numpy_error_curve(D.exponential, D.U_POINT, D.U_EXACT_SLOPE, D.U_WIDTHS, central_difference)
print(f"     dtype               {array_errors.dtype}")
print(f"     shape               {array_errors.shape}")
print(f"     argmin              {int(array_errors.argmin())}")
print(f"     h at that index     {D.U_WIDTHS[int(array_errors.argmin())]:.3e}")
print(f"     same as best_step   {D.U_WIDTHS[int(array_errors.argmin())] == best_central_h}")
assert D.U_WIDTHS[int(array_errors.argmin())] == best_central_h
assert array_errors.shape == (len(D.U_WIDTHS),)
print()
print("   Same numbers, different container. The array form is what makes a")
print("   log-log plot one call, and the U is unmistakable when it is drawn.")
print()

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

"""A zero derivative says the ground is level. It does not say where you are.

Run from inside `examples/`:

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

from __future__ import annotations

import dataset as D
from derivatives import central_difference, classify_stationary_point, second_difference

print("Day 108 / 06 — flat points, and telling them apart")
print()

H = D.STATIONARY_WIDTH
TOL = D.STATIONARY_TOL

# --------------------------------------------------------------------------
print("1. Three places where the slope is zero")
# --------------------------------------------------------------------------
print()
print("     f(x) = (x - 2)**2 + 1   at x = 2     the bottom of a valley")
print("     f(x) = x**3 - 3x        at x = -1    the top of a hill")
print("     f(x) = x**3             at x = 0     neither: a flat step")
print()
print(f"   Measured with the central difference at h = {H:.0e}:")
print()
print("     function            x       f(x)        f'(x) measured")
cases = [
    ("(x - 2)**2 + 1", D.parabola, 2.0),
    ("x**3 - 3x", D.cubic, -1.0),
    ("x**3 - 3x", D.cubic, 1.0),
    ("x**3", D.plain_cube, 0.0),
]
for label, f, x in cases:
    slope = central_difference(f, x, H)
    print(f"     {label:<19} {x:<7.1f} {f(x):<11.4f} {slope:.3e}")
    assert abs(slope) < TOL, (label, x, slope)
print()
print("   Four flat points, four slopes indistinguishable from zero. The first")
print("   derivative has now told you everything it knows, and it has not told")
print("   you which of these is a minimum. That is not a limitation of the")
print("   measurement. It is a limitation of the question.")
print()

# --------------------------------------------------------------------------
print("2. Look at the neighbourhood and the difference is obvious")
# --------------------------------------------------------------------------
print()
print("   Step a little either way from each point and compare the values:")
print()
print("     function            x        f(x - 0.1)   f(x)        f(x + 0.1)   verdict")
for label, f, x in cases:
    left, here, right = f(x - 0.1), f(x), f(x + 0.1)
    if left > here < right:
        verdict = "both sides higher -> minimum"
    elif left < here > right:
        verdict = "both sides lower  -> maximum"
    else:
        verdict = "one of each       -> neither"
    print(f"     {label:<19} {x:<8.1f} {left:<12.5f} {here:<11.5f} {right:<12.5f} {verdict}")
print()
print("   That works, and it is what your eye does when it looks at a graph.")
print("   It is also not a formula, and it needs you to pick 0.1 out of the")
print("   air. The second derivative is the same idea made precise.")
print()

# --------------------------------------------------------------------------
print("3. The second derivative: the rate of change of the rate of change")
# --------------------------------------------------------------------------
print()
print("   The derivative of a function is a function, so it has a derivative")
print("   of its own. Written f''(x), or d2y/dx2. It answers: is the slope")
print("   itself increasing or decreasing as you move right?")
print()
print("   Numerically, taking a central difference of central differences and")
print("   letting the algebra collapse gives one formula:")
print()
print("     f''(x)  ~  ( f(x + h) - 2*f(x) + f(x - h) ) / h**2")
print()
print("   Read it as: how much does the middle sag below the average of its")
print("   two neighbours? Sagging down is positive curvature, a bowl. Bulging")
print("   up is negative curvature, a dome.")
print()
print("     function            x        f'(x)         f''(x)        shape")
for label, f, x in cases:
    slope = central_difference(f, x, H)
    curve = second_difference(f, x, H)
    shape = "bowl" if curve > TOL else ("dome" if curve < -TOL else "flat both ways")
    print(f"     {label:<19} {x:<8.1f} {slope:<13.3e} {curve:<13.6f} {shape}")
print()

parabola_second = second_difference(D.parabola, 2.0, H)
cubic_min_second = second_difference(D.cubic, 1.0, H)
cubic_max_second = second_difference(D.cubic, -1.0, H)
cube_second = second_difference(D.plain_cube, 0.0, H)

assert abs(parabola_second - 2.0) < D.SECOND_TOL, parabola_second
assert abs(cubic_min_second - 6.0) < D.SECOND_TOL, cubic_min_second
assert abs(cubic_max_second - (-6.0)) < D.SECOND_TOL, cubic_max_second
assert abs(cube_second) < D.SECOND_TOL, cube_second
assert cubic_min_second > 0.0 > cubic_max_second

print("   The exact values are 2, 6, -6 and 0, and all four measurements match")
print(f"   them to inside {D.SECOND_TOL:.0e}. Note that the first derivative gave the")
print("   same answer at all four points and the second derivative gave four")
print("   different ones. That is the whole point of computing it.")
print()

# --------------------------------------------------------------------------
print("4. The classification, and the case it refuses to decide")
# --------------------------------------------------------------------------
print()
print("     function            x        classification")
for label, f, x in cases + [("x**3 - 3x", D.cubic, 0.0)]:
    verdict = classify_stationary_point(f, x, H, TOL)
    print(f"     {label:<19} {x:<8.1f} {verdict}")
print()
assert classify_stationary_point(D.parabola, 2.0, H, TOL) == "minimum"
assert classify_stationary_point(D.cubic, 1.0, H, TOL) == "minimum"
assert classify_stationary_point(D.cubic, -1.0, H, TOL) == "maximum"
assert classify_stationary_point(D.plain_cube, 0.0, H, TOL) == "undecided"
assert classify_stationary_point(D.cubic, 0.0, H, TOL) == "not stationary"
print("   'undecided' is not a bug and it is not a failure of the numerics.")
print("   x**3 at 0 is flat and has zero curvature, and so does x**4 at 0.")
print("   The first is a step in a rising slope and the second is a genuine")
print("   minimum, and no amount of second-derivative information separates")
print("   them. A function that reported 'minimum' there would be lying with")
print("   confidence, which is worse than saying it does not know.")
print()
fourth_slope = central_difference(lambda x: x**4, 0.0, H)
fourth_curve = second_difference(lambda x: x**4, 0.0, H)
print(f"     x**4 at 0:  f' = {fourth_slope:.3e}   f'' = {fourth_curve:.3e}   (a real minimum)")
print(f"     x**3 at 0:  f' = {central_difference(D.plain_cube, 0.0, H):.3e}   "
      f"f'' = {cube_second:.3e}   (not a minimum)")
assert abs(fourth_slope) < TOL and abs(fourth_curve) < D.SECOND_TOL
print()
print("   Both readings are zero as far as this method can see -- 2e-8 and 0")
print("   are the same number to a rule whose own rounding noise here is")
print(f"   around {D.SECOND_TOL:.0e}. Same readings, different answers. This is the honest")
print("   boundary of what a second derivative can do.")
print()

# --------------------------------------------------------------------------
print("5. Why any of this matters for training a model")
# --------------------------------------------------------------------------
print()
print("   Training searches for the minimum of a loss function, and the")
print("   derivative is how it knows which way is downhill:")
print()
print("     f'(x) > 0   the function rises to the right, so step LEFT")
print("     f'(x) < 0   the function falls to the right, so step RIGHT")
print("     f'(x) = 0   flat: nothing to learn from the first derivative")
print()
print("   Watch that read out along the parabola, whose minimum is at x = 2:")
print()
print("     x        f(x)        f'(x)        which way is downhill")
for x in [-1.0, 0.5, 1.5, 2.0, 2.5, 4.0]:
    slope = central_difference(D.parabola, x, H)
    if slope > TOL:
        direction = "left"
    elif slope < -TOL:
        direction = "right"
    else:
        direction = "already flat"
    print(f"     {x:<8.1f} {D.parabola(x):<11.4f} {slope:<12.4f} {direction}")
    assert abs(slope - D.parabola_derivative(x)) < D.SECOND_TOL
print()
print("   Every one of those arrows points towards x = 2, and none of them was")
print("   told where x = 2 is. That is gradient descent in one dimension, and")
print("   Day 111 will write the four-line loop that follows the arrows.")
print()

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

"""A corner has no slope -- and the central difference will hand you one anyway.

Run from inside `examples/`:

    ../.venv/bin/python3 07_where_the_derivative_fails.py

This is the script to remember. Everything before it showed a method working.
This one shows the method being confidently, quietly wrong, on the exact shape
you will meet inside every neural network you ever train.
"""

from __future__ import annotations

import dataset as D
from derivatives import (
    backward_difference,
    central_difference,
    forward_difference,
    second_difference,
)

print("Day 108 / 07 — where no derivative exists")
print()

H = D.CORNER_WIDTH

# --------------------------------------------------------------------------
print("1. f(x) = |x| near zero")
# --------------------------------------------------------------------------
print()
print("     x        |x|")
for x in [-0.2, -0.1, 0.0, 0.1, 0.2]:
    print(f"     {x:<8.1f} {D.absolute(x):.1f}")
print()
print("   Two straight lines meeting at a point. To the left of zero the slope")
print("   is -1 everywhere. To the right it is +1 everywhere. At zero it is")
print("   neither, and there is no third answer hiding between them.")
print()
print("   The definition of the derivative asks for a single number that the")
print("   difference quotient settles on as h shrinks -- from BOTH sides. Here")
print("   the two sides settle on different numbers, so the limit does not")
print("   exist, and neither does the derivative. |x| is continuous at zero")
print("   and not differentiable at zero; those are different questions.")
print()

# --------------------------------------------------------------------------
print("2. Ask the three rules anyway")
# --------------------------------------------------------------------------
print()
forward = forward_difference(D.absolute, 0.0, H)
backward = backward_difference(D.absolute, 0.0, H)
central = central_difference(D.absolute, 0.0, H)
print(f"     forward_difference(abs, 0, {H:.0e})    {forward!r}")
print(f"     backward_difference(abs, 0, {H:.0e})   {backward!r}")
print(f"     central_difference(abs, 0, {H:.0e})    {central!r}")
print()
assert forward == D.ABS_FORWARD_AT_ZERO
assert backward == D.ABS_BACKWARD_AT_ZERO
assert central == D.ABS_CENTRAL_AT_ZERO
print("   The two one-sided rules disagree, which is the truth: they are")
print("   reporting the two different slopes that meet here, and their")
print("   disagreement is exactly the reason there is no derivative.")
print()
print("   The central rule returns 0.0. Not an error, not a warning, not a nan.")
print("   Zero, which is the average of -1 and +1, and which is the answer to a")
print("   question nobody asked. It is the average of two slopes rather than")
print("   the slope of anything.")
print()
print("   Worse, 0.0 is a plausible-looking answer. It is what you would get at")
print("   the bottom of a valley, and here it means the opposite -- the")
print("   function is changing as fast as it possibly can in both directions.")
print()
print("   Shrinking h does not help, because nothing is converging:")
print()
print("     h            forward     backward    central")
for h in [1e-2, 1e-5, 1e-8, 1e-11]:
    print(f"     {h:<12.0e} {forward_difference(D.absolute, 0.0, h):<11.1f} "
          f"{backward_difference(D.absolute, 0.0, h):<11.1f} "
          f"{central_difference(D.absolute, 0.0, h):.1f}")
    assert central_difference(D.absolute, 0.0, h) == 0.0
print()
print("   Every row is the same. There is no h small enough to reveal a limit")
print("   that is not there. Compare that with script 02, where the numbers")
print("   visibly settled -- settling is the evidence, and here there is none.")
print()

# --------------------------------------------------------------------------
print("3. The second derivative at a corner is worse still")
# --------------------------------------------------------------------------
print()
print("     h            second_difference(abs, 0, h)")
for h in [1e-2, 1e-3, 1e-5]:
    curve = second_difference(D.absolute, 0.0, h)
    print(f"     {h:<12.0e} {curve:,.1f}")
    assert abs(curve - 2.0 / h) < 1e-6 * (2.0 / h)
print()
print("   It is 2/h, so it grows without limit as h shrinks. A number that")
print("   doubles every time you halve h is not converging on anything, and")
print("   that divergence is a far better warning sign than the first")
print("   derivative's calm 0.0. If a curvature estimate explodes when you")
print("   shrink the step, you are standing on a corner.")
print()

# --------------------------------------------------------------------------
print("4. Why this is not a curiosity: ReLU")
# --------------------------------------------------------------------------
print()
print("   ReLU is max(x, 0): the most widely used activation function in deep")
print("   learning, and Day 102 already met it as a transformation. Its graph")
print("   is flat to the left of zero and a 45-degree line to the right --")
print("   the same corner as |x|, with one arm flattened.")
print()
print("     x        relu(x)")
for x in [-0.2, -0.1, 0.0, 0.1, 0.2]:
    print(f"     {x:<8.1f} {D.relu(x):.1f}")
print()
relu_forward = forward_difference(D.relu, 0.0, H)
relu_backward = backward_difference(D.relu, 0.0, H)
relu_central = central_difference(D.relu, 0.0, H)
print(f"     forward   {relu_forward!r}     the slope on the right")
print(f"     backward  {relu_backward!r}     the slope on the left")
print(f"     central   {relu_central!r}     the average of two slopes that disagree")
assert relu_forward == D.RELU_FORWARD_AT_ZERO
assert relu_backward == D.RELU_BACKWARD_AT_ZERO
assert relu_central == D.RELU_CENTRAL_AT_ZERO
print()
print("   Training a network needs a derivative of ReLU at every input,")
print("   including exactly zero. There is no derivative there, so a framework")
print("   has to choose one of 0 and 1 by convention and carry on. That is an")
print("   engineering decision rather than a mathematical result, and it is")
print("   defensible: an input that is exactly 0.0 in float64 is vanishingly")
print("   rare, and both candidate answers are finite and small.")
print()
print("   This lab does not have a deep-learning framework installed and so")
print("   makes no claim about which value any particular one picks -- check")
print("   your framework's own documentation rather than a course's memory of")
print("   it. What this lab CAN show you is that 0.5, the number the central")
print("   difference produced, is not either of the two defensible choices. If")
print("   you ever use a numerical derivative to check a framework's gradients")
print("   -- which is a real and useful technique -- it will disagree with the")
print("   framework at exactly this point, and the framework will not be wrong.")
print()

# --------------------------------------------------------------------------
print("5. Away from the corner, everything is fine")
# --------------------------------------------------------------------------
print()
print("     x        relu'(x) measured      exact   error")
for x in [-1.0, -0.5, 0.5, 1.0]:
    measured = central_difference(D.relu, x, H)
    exact = 0.0 if x < 0 else 1.0
    print(f"     {x:<8.1f} {measured!r:<22} {exact:<7.1f} {abs(measured - exact):.1e}")
    assert abs(measured - exact) < D.CENTRAL_TOL
print()
print("   On the flat left arm the answer is exactly 0.0, because both sampled")
print("   values are exactly 0.0 and their difference is exactly zero. On the")
print("   sloping right arm it is 1 to within about 2e-12, which is rounding")
print("   error and nothing else: the exact answer 1 is not exactly")
print("   representable as a difference of two nearby float64 values.")
print()
print("   The failure is at one point, not everywhere, and that is why ReLU is")
print("   usable at all. The rule to carry away: a numerical derivative always")
print("   returns a number, and returning a number is not the same as there")
print("   being one. Corners, jumps and vertical tangents all produce")
print("   confident nonsense, and none of them raises an exception.")
print()
print("   The cheapest check costs one extra call: compute the forward and")
print("   backward differences too, and if they disagree by more than the")
print("   tolerance you expect, do not trust the central one.")
print()
print("     function   forward   backward   disagree?   trust the central value?")
for label, f, x in [("x**2", D.square, 3.0), ("|x|", D.absolute, 0.0), ("relu", D.relu, 0.0)]:
    fwd = forward_difference(f, x, H)
    bwd = backward_difference(f, x, H)
    disagree = abs(fwd - bwd) > 1e-3
    print(f"     {label:<10} {fwd:<9.4f} {bwd:<10.4f} {str(disagree):<11} {'no' if disagree else 'yes'}")
print()
assert abs(forward_difference(D.square, 3.0, H) - backward_difference(D.square, 3.0, H)) < 1e-3
assert abs(forward_difference(D.absolute, 0.0, H) - backward_difference(D.absolute, 0.0, H)) > 1e-3
assert abs(forward_difference(D.relu, 0.0, H) - backward_difference(D.relu, 0.0, H)) > 1e-3
print("   That check would have caught both corners here, and it is the same")
print("   two function values the central rule already computed.")
print()

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

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

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

"""The invented data, the functions, the step sizes and the tolerances.

Everything in this lab is computed from the definitions below. Nothing is read
from disk, nothing is downloaded, and no number here was chosen to make a test
pass -- every tolerance is derived in `TOLERANCES` from the two error terms that
actually govern a difference quotient, and the derivation is written out beside
it so you can check the arithmetic yourself.

Read this file. Do not change it: the reference tests compare captured values
against the constants here, so editing one moves the goalposts rather than
fixing anything.
"""

from __future__ import annotations

import math

# ---------------------------------------------------------------------------
# The car, which is where the day starts
# ---------------------------------------------------------------------------

# Invented. A car's distance from a marker post, in metres, sampled once a
# second for six seconds. The numbers were chosen so the arithmetic is doable in
# your head: they are 4 * t**2, so the car is accelerating steadily.
CAR_TIMES_S = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
CAR_DISTANCE_M = [0.0, 4.0, 16.0, 36.0, 64.0, 100.0, 144.0]

# Average speed over the whole six seconds: (144 - 0) / (6 - 0).
CAR_AVERAGE_SPEED_WHOLE_TRIP = 24.0
# Average speed over the fourth second only: (64 - 36) / (4 - 3).
CAR_AVERAGE_SPEED_SECOND_FOUR = 28.0
# The speedometer reading at t = 3 exactly, which is d/dt of 4t**2 = 8t.
CAR_INSTANT_SPEED_AT_3 = 24.0


# ---------------------------------------------------------------------------
# The functions the lab differentiates
# ---------------------------------------------------------------------------


def square(x: float) -> float:
    """f(x) = x**2. Exact derivative 2x. The whole lab's first example."""
    return x * x


def square_derivative(x: float) -> float:
    return 2.0 * x


def cubic(x: float) -> float:
    """f(x) = x**3 - 3x. Stationary at x = -1 (a maximum) and x = +1 (a minimum)."""
    return x * x * x - 3.0 * x


def cubic_derivative(x: float) -> float:
    return 3.0 * x * x - 3.0


def cubic_second_derivative(x: float) -> float:
    return 6.0 * x


def parabola(x: float) -> float:
    """f(x) = (x - 2)**2 + 1. Vertex at x = 2, where the slope is exactly zero."""
    return (x - 2.0) ** 2 + 1.0


def parabola_derivative(x: float) -> float:
    return 2.0 * (x - 2.0)


def plain_cube(x: float) -> float:
    """f(x) = x**3. Flat at x = 0 and yet neither a maximum nor a minimum."""
    return x * x * x


def plain_cube_derivative(x: float) -> float:
    return 3.0 * x * x


def exponential(x: float) -> float:
    """f(x) = e**x. Its own derivative, which is what makes e special."""
    return math.exp(x)


def natural_log(x: float) -> float:
    """f(x) = ln(x). Derivative 1/x."""
    return math.log(x)


def absolute(x: float) -> float:
    """f(x) = |x|. No derivative at 0 -- and the corner ReLU is built from."""
    return abs(x)


def relu(x: float) -> float:
    """max(x, 0). Day 102 met this as a transformation; here it is the corner."""
    return x if x > 0.0 else 0.0


# ---------------------------------------------------------------------------
# The points, the widths and the exact answers
# ---------------------------------------------------------------------------

# Where the shrinking-interval demonstration happens.
SETTLE_POINT = 3.0
SETTLE_EXACT_SLOPE = 6.0  # d/dx of x**2 at x = 3
SETTLE_WIDTHS = [1.0, 0.1, 0.01, 0.001]
# For f(x) = x**2 the average rate over [a, a + h] is exactly 2a + h, so at
# a = 3 the sequence below is 6 + h and can be written down without a computer.
SETTLE_EXPECTED_SLOPES = [7.0, 6.1, 6.01, 6.001]

# Where the U-shaped error curve is measured. e**x at x = 1 is chosen because
# every one of its derivatives is e, which makes both error terms easy to state.
U_POINT = 1.0
U_EXACT_SLOPE = math.e  # 2.718281828459045

# 27 step sizes, one per decade and two thirds, spanning 1e-1 down to 1e-14.
U_WIDTHS = [10.0 ** (-1.0 - 0.5 * k) for k in range(27)]

# The point at which forward and central are compared head to head.
COMPARE_WIDTH = 1e-5

# The step used for the stationary-point work. Small enough that the truncation
# term is invisible, large enough that the second difference's h**2 divisor has
# not started amplifying rounding error.
STATIONARY_WIDTH = 1e-4

# Machine epsilon for float64: the gap between 1.0 and the next float up.
EPSILON = 2.220446049250313e-16


# ---------------------------------------------------------------------------
# The tolerances, and where each one comes from
# ---------------------------------------------------------------------------

# A difference quotient carries two errors that pull in opposite directions.
#
#   TRUNCATION comes from the mathematics: the formula is only the limit's
#   approximation at a finite h. Taylor gives, for f = e**x at x = 1 where every
#   derivative equals e:
#       forward:  |error| ~ (h / 2)  * e
#       central:  |error| ~ (h**2 / 6) * e
#
#   ROUNDING comes from the arithmetic: f(x + h) and f(x - h) are each stored to
#   about EPSILON relative precision, their difference cancels most of their
#   digits, and dividing by h magnifies what is left:
#       either rule: |error| ~ EPSILON * e / h
#
# Add the two, put the numbers in, and every tolerance below follows. None was
# reached by running a test and enlarging the number until it went green.

# Central difference of e**x at x = 1 with h = 1e-5:
#   truncation ~ (1e-10 / 6) * e   = 4.5e-11
#   rounding   ~ 2.22e-16 * e / 1e-5 = 6.0e-11
#   sum        ~ 1.1e-10            -> allow 1e-9, roughly nine times the bound
CENTRAL_TOL = 1e-9

# Forward difference of e**x at x = 1 with h = 1e-5:
#   truncation ~ (1e-5 / 2) * e = 1.36e-5
#   rounding   ~ 6.0e-11, negligible beside it
#   sum        ~ 1.4e-5         -> allow 1e-4, seven times the bound
FORWARD_TOL = 1e-4

# The average-rate sequence over [3, 3 + h] for f = x**2. Each term is 2a + h
# computed in float64 from numbers of order 10, so only a few units in the last
# place of about 1e-15 are available to go wrong.
EXACT_TOL = 1e-12

# The second difference of a cubic is exact in the mathematics -- the h**2 term
# in its Taylor expansion is the answer and there is no h**4 term to leave
# behind -- so only rounding is in play. With h = 1e-4 and |f| of order 2:
#   rounding ~ 4 * EPSILON * |f| / h**2 = 4 * 2.22e-16 * 2 / 1e-8 = 1.8e-7
#   -> allow 1e-5, about fifty times the bound
SECOND_TOL = 1e-5

# A stationary point found with the central difference at h = 1e-4 on the cubic:
# the truncation term is exactly h**2 * f'''/6 = h**2, which is 1e-8.
#   -> allow 1e-6, a hundred times the bound
STATIONARY_TOL = 1e-6


# ---------------------------------------------------------------------------
# The derivative rules, stated as facts and checked numerically in the lab
# ---------------------------------------------------------------------------

# (name, f, exact f', a point to check it at)
RULE_CASES = [
    ("constant: d/dx of 7 is 0", lambda x: 7.0, lambda x: 0.0, 2.0),
    ("power: d/dx of x**2 is 2x", square, square_derivative, 3.0),
    ("power: d/dx of x**5 is 5x**4", lambda x: x**5, lambda x: 5.0 * x**4, 1.5),
    ("power: d/dx of 1/x is -1/x**2", lambda x: 1.0 / x, lambda x: -1.0 / (x * x), 2.0),
    ("constant multiple: d/dx of 5x**2 is 10x", lambda x: 5.0 * x * x, lambda x: 10.0 * x, 3.0),
    ("sum: d/dx of x**2 + x**3 is 2x + 3x**2", lambda x: x**2 + x**3, lambda x: 2.0 * x + 3.0 * x**2, 2.0),
    ("exponential: d/dx of e**x is e**x", exponential, exponential, 1.0),
    ("logarithm: d/dx of ln(x) is 1/x", natural_log, lambda x: 1.0 / x, 4.0),
]

# The exact slopes those eight cases must produce, written out so a reader can
# check them by hand rather than by rerunning the lab.
RULE_EXPECTED = [0.0, 6.0, 25.3125, -0.25, 30.0, 16.0, math.e, 0.25]

# The eight rule cases are checked with the central difference at h = 1e-5. The
# worst truncation term among them belongs to x**5 at 1.5, where the third
# derivative is 60 * 1.5**2 = 135:
#   truncation ~ (h**2 / 6) * 135 = (1e-10 / 6) * 135 = 2.25e-9
#   -> allow 1e-8, about four times the bound
RULE_TOL = 1e-8

# The corner cases at x = 0, where no derivative exists.
CORNER_WIDTH = 1e-5
# |x|: the one-sided slopes are -1 and +1 and disagree, which is the whole
# reason there is no derivative. The central difference averages them to zero
# and reports that zero with total confidence.
ABS_FORWARD_AT_ZERO = 1.0
ABS_BACKWARD_AT_ZERO = -1.0
ABS_CENTRAL_AT_ZERO = 0.0
# max(x, 0): the one-sided slopes are 0 and 1, so the central difference gives
# their average, 0.5. Deep-learning frameworks do not use 0.5; they pick one of
# the one-sided values and move on, and that choice is a convention rather than
# a theorem.
RELU_FORWARD_AT_ZERO = 1.0
RELU_BACKWARD_AT_ZERO = 0.0
RELU_CENTRAL_AT_ZERO = 0.5

"""The reference implementation: rates of change, computed from first principles.

Ten functions, all built on one idea -- a rate is a difference divided by the
interval it happened over -- and one NumPy-based helper at the end so the
from-scratch version can be checked against the library.

Read `starter/derivatives.py` and write your own before you read this one.
"""

from __future__ import annotations

from collections.abc import Callable, Sequence

import numpy as np

Function = Callable[[float], float]


# ---------------------------------------------------------------------------
# 1. Average rate of change: rise over run, with real numbers
# ---------------------------------------------------------------------------


def average_rate(f: Function, a: float, b: float) -> float:
    """The average rate of change of f between a and b.

    Rise over run, and nothing more: how much the output changed, divided by
    how much the input changed. If f is distance in metres and the inputs are
    seconds, this is metres per second -- an average speed over the interval,
    exactly what a stopwatch and a tape measure would give you.

    Raises ZeroDivisionError when a == b, which is the honest thing to do: the
    question "how fast, over an interval of no width" has no answer, and the
    whole of today is about approaching that question rather than asking it.
    """
    if a == b:
        raise ZeroDivisionError(
            "average_rate needs an interval with width; a and b are both "
            f"{a!r}. The rate 'right here' is what the derivative is for."
        )
    return (f(b) - f(a)) / (b - a)


def shrinking_slopes(f: Function, a: float, widths: Sequence[float]) -> list[float]:
    """The average rate over [a, a + h] for each h, in the order given.

    Feed it widths that get smaller and watch the returned numbers settle. That
    settling is the limit, met as an observation rather than as a definition.
    """
    return [average_rate(f, a, a + h) for h in widths]


# ---------------------------------------------------------------------------
# 2. The two difference quotients
# ---------------------------------------------------------------------------


def forward_difference(f: Function, x: float, h: float) -> float:
    """Slope of the secant from x to x + h. The definition, stopped early.

    Error shrinks like h: halve the step and you roughly halve the error.
    """
    return (f(x + h) - f(x)) / h


def backward_difference(f: Function, x: float, h: float) -> float:
    """Slope of the secant from x - h to x. The forward rule facing the other way."""
    return (f(x) - f(x - h)) / h


def central_difference(f: Function, x: float, h: float) -> float:
    """Slope of the secant from x - h to x + h, straddling the point.

    This is the average of the forward and backward differences, and the
    averaging is why it is so much better: the two rules lean off the tangent in
    opposite directions by almost exactly the same amount, so their leading
    errors cancel. Error shrinks like h**2 -- halve the step and the error
    quarters -- for one extra function call over the forward rule, and none at
    all over computing forward and backward separately.
    """
    return (f(x + h) - f(x - h)) / (2.0 * h)


def second_difference(f: Function, x: float, h: float) -> float:
    """The rate of change of the rate of change: an approximate f''(x).

    Built by taking a central difference of central differences and letting the
    algebra collapse. Its sign is curvature: positive is a bowl, negative is a
    dome, and that is what tells a minimum from a maximum when the first
    derivative has said only 'flat'.

    Note the h**2 in the divisor. It magnifies rounding error far harder than
    the first-difference rules do, so the useful range of h is both narrower and
    larger here.
    """
    return (f(x + h) - 2.0 * f(x) + f(x - h)) / (h * h)


# ---------------------------------------------------------------------------
# 3. Measuring how wrong the approximation is
# ---------------------------------------------------------------------------


def error_curve(
    f: Function,
    x: float,
    exact_slope: float,
    widths: Sequence[float],
    rule: Callable[[Function, float, float], float],
) -> list[float]:
    """Absolute error of `rule` against a known exact slope, one entry per width.

    Only usable when you already know the right answer, which is exactly why the
    lab measures it on e**x rather than on something interesting: the point is
    to see the shape of the error, and you cannot see the shape of an error you
    cannot compute.
    """
    return [abs(rule(f, x, h) - exact_slope) for h in widths]


def best_step(widths: Sequence[float], errors: Sequence[float]) -> tuple[float, float]:
    """The (width, error) pair with the smallest error. The bottom of the U.

    Ties go to the first, which for a descending list of widths means the
    largest h that achieves the minimum -- the conservative choice, since a
    larger step sits further from the cancellation cliff.
    """
    if len(widths) != len(errors):
        raise ValueError(
            f"widths and errors must be the same length; got {len(widths)} and {len(errors)}"
        )
    if not widths:
        raise ValueError("best_step needs at least one width")
    index = min(range(len(errors)), key=lambda i: errors[i])
    return (widths[index], errors[index])


def is_u_shaped(errors: Sequence[float]) -> bool:
    """True when the errors fall to a single minimum and then rise again.

    Deliberately tolerant of the wobble at the bottom: rounding error is a
    random walk, not a smooth curve, so this asks only that the error at the
    large-h end and the error at the small-h end are both meaningfully worse
    than the best one in the middle, and that the minimum is not at either end.
    """
    if len(errors) < 3:
        return False
    index = min(range(len(errors)), key=lambda i: errors[i])
    if index == 0 or index == len(errors) - 1:
        return False
    best = errors[index]
    return errors[0] > 10.0 * best and errors[-1] > 10.0 * best


# ---------------------------------------------------------------------------
# 4. What a zero derivative does and does not tell you
# ---------------------------------------------------------------------------


def tangent_at(f: Function, x: float, h: float) -> tuple[float, float]:
    """(slope, intercept) of the tangent line to f at x, estimated centrally.

    The tangent is the line through (x, f(x)) with the derivative as its slope,
    so the intercept follows from y = mx + c rearranged: c = f(x) - m*x.
    """
    slope = central_difference(f, x, h)
    return (slope, f(x) - slope * x)


def classify_stationary_point(f: Function, x: float, h: float, tol: float) -> str:
    """Name what kind of point x is, from the first and second derivatives.

    Returns one of:
      'not stationary' -- the first derivative is not zero within tol
      'minimum'        -- flat, and curving upward
      'maximum'        -- flat, and curving downward
      'undecided'      -- flat, and the second derivative is zero too

    That last case is the honest one and the reason this function exists. A zero
    first derivative says the ground is level; it does not say whether you are
    at the bottom of a valley, the top of a hill, or on a flat step partway down
    a slope. The second derivative resolves two of those three, and when it is
    also zero it resolves nothing -- x**3 at 0 and x**4 at 0 are both flat with
    zero curvature and are a step and a minimum respectively.
    """
    first = central_difference(f, x, h)
    if abs(first) > tol:
        return "not stationary"
    second = second_difference(f, x, h)
    if second > tol:
        return "minimum"
    if second < -tol:
        return "maximum"
    return "undecided"


# ---------------------------------------------------------------------------
# 5. The same job, handed to NumPy
# ---------------------------------------------------------------------------


def numpy_gradient_slope(f: Function, x: float, h: float) -> float:
    """The derivative at x via numpy.gradient over a three-point sample.

    `np.gradient` differentiates SAMPLES rather than a function: you hand it
    values you already have and it returns a derivative estimate at every one of
    them. On the interior points it uses the central difference, which is why
    the middle of a three-point sample straddling x agrees with
    `central_difference` to the last bit. On the two ends it has nothing on one
    side, so it falls back to a one-sided rule -- which is exactly the forward
    or backward difference, and exactly as much worse.

    The spacing is passed as the scalar `h`. Handing `np.gradient` an array of
    coordinates instead is algebraically the same request and takes a different
    route through the arithmetic, so it lands a few units in the last place away;
    `numpy_gradient_slope_from_coordinates` below is that version, kept so the
    difference can be measured rather than argued about.
    """
    ys = np.array([f(x - h), f(x), f(x + h)], dtype=np.float64)
    return float(np.gradient(ys, h)[1])


def numpy_gradient_slope_from_coordinates(f: Function, x: float, h: float) -> float:
    """The same three-point estimate, with coordinates passed instead of spacing.

    Kept only to show that "the same formula" is a claim about the mathematics.
    With unevenly spaced coordinates NumPy must use a general weighted rule, and
    it uses that rule even when the coordinates happen to be evenly spaced -- so
    this returns a number a few units in the last place from the one above.
    """
    xs = np.array([x - h, x, x + h], dtype=np.float64)
    ys = np.array([f(v) for v in xs], dtype=np.float64)
    return float(np.gradient(ys, xs)[1])


def numpy_error_curve(
    f: Function,
    x: float,
    exact_slope: float,
    widths: Sequence[float],
    rule: Callable[[Function, float, float], float],
) -> np.ndarray:
    """`error_curve` as a float64 array, for plotting and for argmin.

    Same numbers, different container. The array form is what makes
    `errors.argmin()` and the log-log plot in the lesson one line each.
    """
    return np.array(error_curve(f, x, exact_slope, widths, rule), dtype=np.float64)

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

Run from the lab directory:

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

Nothing here reads source code or checks that a function exists. Every test
calls something and compares the result against a number that was either
derived by hand or derived from the error analysis written out in dataset.py.
"""

from __future__ import annotations

import math

import numpy as np
import pytest

import dataset as D
from derivatives import (
    average_rate,
    backward_difference,
    best_step,
    central_difference,
    classify_stationary_point,
    error_curve,
    forward_difference,
    is_u_shaped,
    numpy_error_curve,
    numpy_gradient_slope,
    numpy_gradient_slope_from_coordinates,
    second_difference,
    shrinking_slopes,
    tangent_at,
)


# ---------------------------------------------------------------------------
# Average rate of change
# ---------------------------------------------------------------------------


def car_distance(t: float) -> float:
    return 4.0 * t * t


def test_car_table_matches_the_formula():
    assert [car_distance(t) for t in D.CAR_TIMES_S] == D.CAR_DISTANCE_M


def test_average_speed_over_the_whole_trip():
    assert average_rate(car_distance, 0.0, 6.0) == D.CAR_AVERAGE_SPEED_WHOLE_TRIP


def test_average_speed_over_the_fourth_second():
    assert average_rate(car_distance, 3.0, 4.0) == D.CAR_AVERAGE_SPEED_SECOND_FOUR


@pytest.mark.parametrize(
    ("a", "b", "expected"),
    [(0.0, 1.0, 4.0), (1.0, 2.0, 12.0), (2.0, 3.0, 20.0), (3.0, 4.0, 28.0), (4.0, 5.0, 36.0), (5.0, 6.0, 44.0)],
)
def test_second_by_second_average_speeds(a, b, expected):
    assert average_rate(car_distance, a, b) == expected


def test_average_rate_is_symmetric_in_its_endpoints():
    assert average_rate(D.square, 1.0, 4.0) == average_rate(D.square, 4.0, 1.0)


def test_average_rate_over_a_zero_width_interval_refuses():
    with pytest.raises(ZeroDivisionError):
        average_rate(car_distance, 3.0, 3.0)


def test_the_refusal_explains_itself():
    with pytest.raises(ZeroDivisionError) as caught:
        average_rate(car_distance, 3.0, 3.0)
    assert "interval with width" in str(caught.value)


def test_average_rate_of_a_straight_line_is_the_same_everywhere():
    line = lambda x: 3.0 * x + 5.0  # noqa: E731 - deliberately inline
    assert average_rate(line, 0.0, 1.0) == 3.0
    assert average_rate(line, -100.0, 100.0) == 3.0


def test_instantaneous_speed_sits_between_its_neighbouring_averages():
    before = average_rate(car_distance, 2.0, 3.0)
    after = average_rate(car_distance, 3.0, 4.0)
    assert before < D.CAR_INSTANT_SPEED_AT_3 < after


# ---------------------------------------------------------------------------
# The shrinking interval
# ---------------------------------------------------------------------------


def test_shrinking_slopes_returns_one_value_per_width():
    slopes = shrinking_slopes(D.square, D.SETTLE_POINT, D.SETTLE_WIDTHS)
    assert len(slopes) == len(D.SETTLE_WIDTHS)


@pytest.mark.parametrize(("index", "expected"), list(enumerate(D.SETTLE_EXPECTED_SLOPES)))
def test_each_shrinking_slope_is_six_plus_h(index, expected):
    slopes = shrinking_slopes(D.square, D.SETTLE_POINT, D.SETTLE_WIDTHS)
    assert abs(slopes[index] - expected) < D.EXACT_TOL


def test_the_shrinking_sequence_converges_on_the_exact_derivative():
    slopes = shrinking_slopes(D.square, D.SETTLE_POINT, D.SETTLE_WIDTHS)
    gaps = [abs(s - D.SETTLE_EXACT_SLOPE) for s in slopes]
    assert gaps == sorted(gaps, reverse=True)
    assert gaps[-1] < 0.002


def test_the_sequence_approaches_from_above_for_a_convex_function():
    slopes = shrinking_slopes(D.square, D.SETTLE_POINT, D.SETTLE_WIDTHS)
    assert all(s > D.SETTLE_EXACT_SLOPE for s in slopes)


def test_shrinking_from_the_left_approaches_the_same_number():
    widths = [-w for w in D.SETTLE_WIDTHS]
    slopes = shrinking_slopes(D.square, D.SETTLE_POINT, widths)
    assert all(s < D.SETTLE_EXACT_SLOPE for s in slopes)
    assert abs(slopes[-1] - D.SETTLE_EXACT_SLOPE) < 0.002


def test_the_car_sequence_settles_on_twenty_four():
    slopes = shrinking_slopes(car_distance, 3.0, [1.0, 0.1, 0.01, 0.001])
    assert abs(slopes[-1] - D.CAR_INSTANT_SPEED_AT_3) < 0.005


def test_tangent_line_at_three_is_six_x_minus_nine():
    slope, intercept = tangent_at(D.square, 3.0, D.COMPARE_WIDTH)
    assert abs(slope - 6.0) < D.CENTRAL_TOL
    assert abs(intercept + 9.0) < 1e-8


def test_the_tangent_line_touches_the_curve_at_the_point():
    slope, intercept = tangent_at(D.square, 3.0, D.COMPARE_WIDTH)
    assert abs((slope * 3.0 + intercept) - D.square(3.0)) < 1e-9


# ---------------------------------------------------------------------------
# The rules
# ---------------------------------------------------------------------------


@pytest.mark.parametrize(
    ("case", "expected"),
    [(case, expected) for case, expected in zip(D.RULE_CASES, D.RULE_EXPECTED)],
    ids=[case[0] for case in D.RULE_CASES],
)
def test_each_rule_states_the_documented_exact_slope(case, expected):
    _, _, exact_derivative, x = case
    assert abs(exact_derivative(x) - expected) < 1e-15


@pytest.mark.parametrize("case", D.RULE_CASES, ids=[case[0] for case in D.RULE_CASES])
def test_each_rule_agrees_with_the_arithmetic(case):
    _, f, exact_derivative, x = case
    measured = central_difference(f, x, D.COMPARE_WIDTH)
    assert abs(measured - exact_derivative(x)) < D.RULE_TOL


def test_the_constant_rule_gives_exactly_zero():
    assert central_difference(lambda x: 7.0, 2.0, D.COMPARE_WIDTH) == 0.0


def test_the_sum_rule_is_the_sum_of_the_parts():
    at = 2.0
    both = central_difference(lambda x: x**2 + x**3, at, D.COMPARE_WIDTH)
    separately = central_difference(D.square, at, D.COMPARE_WIDTH) + central_difference(
        D.plain_cube, at, D.COMPARE_WIDTH
    )
    assert abs(both - separately) < D.RULE_TOL


def test_the_constant_multiple_rule_scales_the_slope():
    at = 3.0
    plain = central_difference(D.square, at, D.COMPARE_WIDTH)
    scaled = central_difference(lambda x: 5.0 * x * x, at, D.COMPARE_WIDTH)
    assert abs(scaled - 5.0 * plain) < D.RULE_TOL


@pytest.mark.parametrize(("base", "expected"), [(2.0, math.log(2.0)), (3.0, math.log(3.0)), (10.0, math.log(10.0))])
def test_the_slope_of_b_to_the_x_at_zero_is_the_natural_log_of_b(base, expected):
    measured = central_difference(lambda x: base**x, 0.0, D.COMPARE_WIDTH)
    assert abs(measured - expected) < 1e-8


def test_e_is_the_base_whose_slope_at_zero_is_one():
    measured = central_difference(D.exponential, 0.0, D.COMPARE_WIDTH)
    assert abs(measured - 1.0) < D.CENTRAL_TOL


@pytest.mark.parametrize("x", [0.0, 1.0, 2.5])
def test_the_exponential_is_proportional_to_its_own_derivative_with_constant_one(x):
    ratio = central_difference(D.exponential, x, D.COMPARE_WIDTH) / D.exponential(x)
    assert abs(ratio - 1.0) < 1e-9


@pytest.mark.parametrize("x", [0.5, 1.0, 4.0, 10.0])
def test_the_derivative_of_ln_is_one_over_x(x):
    measured = central_difference(D.natural_log, x, D.COMPARE_WIDTH)
    assert abs(measured - 1.0 / x) < 1e-7


# ---------------------------------------------------------------------------
# Forward, backward, central
# ---------------------------------------------------------------------------


@pytest.mark.parametrize("h", [1.0, 0.1, 0.01, 0.001])
def test_forward_difference_of_x_squared_is_exactly_six_plus_h(h):
    assert abs(forward_difference(D.square, 3.0, h) - (6.0 + h)) < 1e-9


@pytest.mark.parametrize("h", [1.0, 0.1, 0.01, 0.001])
def test_backward_difference_of_x_squared_is_exactly_six_minus_h(h):
    assert abs(backward_difference(D.square, 3.0, h) - (6.0 - h)) < 1e-9


@pytest.mark.parametrize("h", [1.0, 0.1, 0.01, 0.001])
def test_central_difference_is_exact_on_a_quadratic(h):
    assert abs(central_difference(D.square, 3.0, h) - 6.0) < 1e-11


def test_central_is_the_average_of_forward_and_backward():
    h = 0.1
    average = (forward_difference(D.exponential, 1.0, h) + backward_difference(D.exponential, 1.0, h)) / 2.0
    assert abs(average - central_difference(D.exponential, 1.0, h)) < 1e-14


def test_forward_difference_meets_its_documented_tolerance():
    error = abs(forward_difference(D.exponential, 1.0, D.COMPARE_WIDTH) - math.e)
    assert error < D.FORWARD_TOL


def test_central_difference_meets_its_documented_tolerance():
    error = abs(central_difference(D.exponential, 1.0, D.COMPARE_WIDTH) - math.e)
    assert error < D.CENTRAL_TOL


def test_central_beats_forward_by_at_least_a_thousandfold_at_the_same_step():
    forward_error = abs(forward_difference(D.exponential, 1.0, D.COMPARE_WIDTH) - math.e)
    central_error = abs(central_difference(D.exponential, 1.0, D.COMPARE_WIDTH) - math.e)
    assert central_error * 1000.0 < forward_error


@pytest.mark.parametrize("h", [1e-1, 1e-2, 1e-3, 1e-4, 1e-5])
def test_central_beats_forward_at_every_sensible_step(h):
    assert abs(central_difference(D.exponential, 1.0, h) - math.e) < abs(
        forward_difference(D.exponential, 1.0, h) - math.e
    )


@pytest.mark.parametrize("h", [1e-1, 1e-2, 1e-3])
def test_forward_error_falls_like_h(h):
    error = abs(forward_difference(D.exponential, 1.0, h) - math.e)
    predicted = (h / 2.0) * math.e
    assert 0.9 < error / predicted < 1.1


@pytest.mark.parametrize("h", [1e-2, 1e-3, 1e-4])
def test_central_error_falls_like_h_squared(h):
    error = abs(central_difference(D.exponential, 1.0, h) - math.e)
    predicted = (h * h / 6.0) * math.e
    assert 0.99 < error / predicted < 1.01


def test_halving_the_step_quarters_the_central_error():
    coarse = abs(central_difference(D.exponential, 1.0, 1e-2) - math.e)
    fine = abs(central_difference(D.exponential, 1.0, 5e-3) - math.e)
    assert 3.9 < coarse / fine < 4.1


def test_halving_the_step_halves_the_forward_error():
    coarse = abs(forward_difference(D.exponential, 1.0, 1e-2) - math.e)
    fine = abs(forward_difference(D.exponential, 1.0, 5e-3) - math.e)
    assert 1.9 < coarse / fine < 2.1


# ---------------------------------------------------------------------------
# The U-shaped error curve
# ---------------------------------------------------------------------------


@pytest.fixture(scope="module")
def forward_errors():
    return error_curve(D.exponential, D.U_POINT, D.U_EXACT_SLOPE, D.U_WIDTHS, forward_difference)


@pytest.fixture(scope="module")
def central_errors():
    return error_curve(D.exponential, D.U_POINT, D.U_EXACT_SLOPE, D.U_WIDTHS, central_difference)


def test_the_grid_spans_the_documented_range():
    assert len(D.U_WIDTHS) == 27
    assert D.U_WIDTHS[0] == 1e-1
    assert D.U_WIDTHS[-1] == 1e-14
    assert D.U_WIDTHS == sorted(D.U_WIDTHS, reverse=True)


def test_error_curve_returns_one_error_per_width(forward_errors):
    assert len(forward_errors) == len(D.U_WIDTHS)


def test_the_forward_error_curve_is_u_shaped(forward_errors):
    assert is_u_shaped(forward_errors)


def test_the_central_error_curve_is_u_shaped(central_errors):
    assert is_u_shaped(central_errors)


def test_the_forward_minimum_is_not_at_either_end(forward_errors):
    index = forward_errors.index(min(forward_errors))
    assert 0 < index < len(forward_errors) - 1


def test_the_central_minimum_is_not_at_either_end(central_errors):
    index = central_errors.index(min(central_errors))
    assert 0 < index < len(central_errors) - 1


def test_the_largest_step_is_far_worse_than_the_best(central_errors):
    best = min(central_errors)
    assert central_errors[0] > 100.0 * best


def test_the_smallest_step_is_far_worse_than_the_best(central_errors):
    best = min(central_errors)
    assert central_errors[-1] > 100.0 * best


def test_the_error_falls_monotonically_down_the_truncation_side(central_errors):
    first_six = central_errors[:6]
    assert first_six == sorted(first_six, reverse=True)


def test_the_error_rises_over_the_last_stretch(central_errors):
    assert central_errors[-1] > central_errors[-6]


def test_best_step_finds_the_bottom(central_errors):
    h, error = best_step(D.U_WIDTHS, central_errors)
    assert error == min(central_errors)
    assert h == D.U_WIDTHS[central_errors.index(error)]


def test_the_best_central_step_is_in_the_documented_band(central_errors):
    h, _ = best_step(D.U_WIDTHS, central_errors)
    assert 1e-7 <= h <= 1e-4


def test_the_best_forward_step_is_in_the_documented_band(forward_errors):
    h, _ = best_step(D.U_WIDTHS, forward_errors)
    assert 1e-9 <= h <= 1e-6


def test_the_best_central_step_beats_the_best_forward_step(forward_errors, central_errors):
    _, forward_best = best_step(D.U_WIDTHS, forward_errors)
    _, central_best = best_step(D.U_WIDTHS, central_errors)
    assert central_best < forward_best


def test_the_measured_optimum_matches_the_balance_prediction(forward_errors, central_errors):
    forward_h, _ = best_step(D.U_WIDTHS, forward_errors)
    central_h, _ = best_step(D.U_WIDTHS, central_errors)
    assert 0.1 < forward_h / math.sqrt(2.0 * D.EPSILON) < 10.0
    assert 0.1 < central_h / (3.0 * D.EPSILON) ** (1.0 / 3.0) < 10.0


def test_a_thousandth_beats_a_millionth_of_a_millionth(central_errors):
    coarse = central_errors[D.U_WIDTHS.index(1e-3)]
    absurd = central_errors[D.U_WIDTHS.index(1e-12)]
    assert absurd > 100.0 * coarse


def test_an_absurdly_small_step_returns_exactly_zero():
    assert forward_difference(D.exponential, 1.0, 1e-300) == 0.0
    assert central_difference(D.exponential, 1.0, 1e-300) == 0.0


def test_the_two_sampled_values_become_the_same_float():
    assert math.exp(1.0 + 1e-300) == math.exp(1.0)


def test_is_u_shaped_rejects_a_monotone_sequence():
    assert not is_u_shaped([100.0, 10.0, 1.0, 0.1])


def test_is_u_shaped_rejects_a_sequence_that_is_too_short():
    assert not is_u_shaped([1.0, 0.5])


def test_is_u_shaped_accepts_a_clear_u():
    assert is_u_shaped([100.0, 1.0, 0.01, 1.0, 100.0])


def test_best_step_rejects_mismatched_lengths():
    with pytest.raises(ValueError):
        best_step([1.0, 2.0], [1.0])


def test_best_step_rejects_an_empty_grid():
    with pytest.raises(ValueError):
        best_step([], [])


def test_the_array_error_curve_holds_the_same_numbers(central_errors):
    array = numpy_error_curve(D.exponential, D.U_POINT, D.U_EXACT_SLOPE, D.U_WIDTHS, central_difference)
    assert array.dtype == np.float64
    assert array.shape == (len(D.U_WIDTHS),)
    assert array.tolist() == central_errors


def test_argmin_agrees_with_best_step(central_errors):
    array = numpy_error_curve(D.exponential, D.U_POINT, D.U_EXACT_SLOPE, D.U_WIDTHS, central_difference)
    h, _ = best_step(D.U_WIDTHS, central_errors)
    assert D.U_WIDTHS[int(array.argmin())] == h


# ---------------------------------------------------------------------------
# Stationary points and curvature
# ---------------------------------------------------------------------------

H = D.STATIONARY_WIDTH
TOL = D.STATIONARY_TOL


@pytest.mark.parametrize(
    ("f", "x"),
    [(D.parabola, 2.0), (D.cubic, 1.0), (D.cubic, -1.0), (D.plain_cube, 0.0)],
    ids=["parabola vertex", "cubic minimum", "cubic maximum", "cubic flat step"],
)
def test_the_first_derivative_is_zero_at_every_stationary_point(f, x):
    assert abs(central_difference(f, x, H)) < TOL


def test_the_parabola_vertex_gives_exactly_zero():
    assert central_difference(D.parabola, 2.0, H) == 0.0


def test_the_cubic_is_not_stationary_at_the_origin():
    assert abs(central_difference(D.cubic, 0.0, H)) > 1.0


def test_the_second_derivative_of_the_parabola_is_two():
    assert abs(second_difference(D.parabola, 2.0, H) - 2.0) < D.SECOND_TOL


def test_the_second_derivative_at_the_cubic_minimum_is_six():
    assert abs(second_difference(D.cubic, 1.0, H) - 6.0) < D.SECOND_TOL


def test_the_second_derivative_at_the_cubic_maximum_is_minus_six():
    assert abs(second_difference(D.cubic, -1.0, H) + 6.0) < D.SECOND_TOL


def test_the_second_derivative_at_the_flat_step_is_zero():
    assert abs(second_difference(D.plain_cube, 0.0, H)) < D.SECOND_TOL


def test_the_second_derivative_separates_the_minimum_from_the_maximum():
    assert second_difference(D.cubic, 1.0, H) > 0.0 > second_difference(D.cubic, -1.0, H)


def test_the_first_derivative_does_not_separate_them():
    at_minimum = central_difference(D.cubic, 1.0, H)
    at_maximum = central_difference(D.cubic, -1.0, H)
    assert abs(at_minimum) < TOL and abs(at_maximum) < TOL
    assert abs(at_minimum - at_maximum) < TOL


def test_the_second_difference_matches_the_exact_second_derivative_along_the_cubic():
    for x in [-2.0, -0.5, 0.5, 2.0]:
        assert abs(second_difference(D.cubic, x, H) - D.cubic_second_derivative(x)) < D.SECOND_TOL


@pytest.mark.parametrize(
    ("f", "x", "expected"),
    [
        (D.parabola, 2.0, "minimum"),
        (D.cubic, 1.0, "minimum"),
        (D.cubic, -1.0, "maximum"),
        (D.plain_cube, 0.0, "undecided"),
        (D.cubic, 0.0, "not stationary"),
        (D.square, 0.0, "minimum"),
        (D.exponential, 0.0, "not stationary"),
    ],
    ids=["parabola", "cubic min", "cubic max", "flat step", "sloping", "x squared", "exponential"],
)
def test_classification(f, x, expected):
    assert classify_stationary_point(f, x, H, TOL) == expected


def test_x_to_the_fourth_at_zero_is_also_undecided():
    assert classify_stationary_point(lambda x: x**4, 0.0, H, TOL) == "undecided"


def test_a_minimum_and_a_step_can_be_indistinguishable():
    step = (central_difference(D.plain_cube, 0.0, H), second_difference(D.plain_cube, 0.0, H))
    minimum = (central_difference(lambda x: x**4, 0.0, H), second_difference(lambda x: x**4, 0.0, H))
    assert abs(step[0] - minimum[0]) < TOL
    assert abs(step[1] - minimum[1]) < D.SECOND_TOL


@pytest.mark.parametrize("x", [-1.0, 0.5, 1.5, 2.5, 4.0])
def test_the_sign_of_the_derivative_points_towards_the_minimum(x):
    slope = central_difference(D.parabola, x, H)
    # A small step against the sign of the slope must land nearer the minimum.
    # The step is 0.1 rather than 1.0 because a step larger than the distance
    # to the minimum overshoots it -- which is Day 111's learning-rate problem,
    # arriving three days early.
    downhill = x - 0.1 * math.copysign(1.0, slope)
    assert abs(downhill - 2.0) < abs(x - 2.0)


@pytest.mark.parametrize("x", [-1.0, 0.5, 1.5, 2.5, 4.0])
def test_the_measured_parabola_slope_matches_the_exact_one(x):
    assert abs(central_difference(D.parabola, x, H) - D.parabola_derivative(x)) < D.SECOND_TOL


# ---------------------------------------------------------------------------
# Where no derivative exists
# ---------------------------------------------------------------------------


def test_the_forward_difference_of_abs_at_zero_is_plus_one():
    assert forward_difference(D.absolute, 0.0, D.CORNER_WIDTH) == D.ABS_FORWARD_AT_ZERO


def test_the_backward_difference_of_abs_at_zero_is_minus_one():
    assert backward_difference(D.absolute, 0.0, D.CORNER_WIDTH) == D.ABS_BACKWARD_AT_ZERO


def test_the_central_difference_of_abs_at_zero_is_zero():
    assert central_difference(D.absolute, 0.0, D.CORNER_WIDTH) == D.ABS_CENTRAL_AT_ZERO


@pytest.mark.parametrize("h", [1e-2, 1e-5, 1e-8, 1e-11, 1e-14])
def test_shrinking_h_never_reveals_a_limit_at_the_corner(h):
    assert central_difference(D.absolute, 0.0, h) == 0.0
    assert forward_difference(D.absolute, 0.0, h) == 1.0
    assert backward_difference(D.absolute, 0.0, h) == -1.0


def test_the_one_sided_rules_disagree_at_the_corner():
    gap = abs(
        forward_difference(D.absolute, 0.0, D.CORNER_WIDTH) - backward_difference(D.absolute, 0.0, D.CORNER_WIDTH)
    )
    assert gap == 2.0


def test_the_one_sided_rules_agree_where_a_derivative_exists():
    gap = abs(forward_difference(D.square, 3.0, D.CORNER_WIDTH) - backward_difference(D.square, 3.0, D.CORNER_WIDTH))
    assert gap < 1e-3


@pytest.mark.parametrize("h", [1e-2, 1e-3, 1e-5])
def test_the_second_difference_at_the_corner_diverges_like_two_over_h(h):
    curve = second_difference(D.absolute, 0.0, h)
    assert abs(curve - 2.0 / h) < 1e-6 * (2.0 / h)


def test_the_corner_curvature_grows_as_the_step_shrinks():
    coarse = second_difference(D.absolute, 0.0, 1e-2)
    fine = second_difference(D.absolute, 0.0, 1e-4)
    assert fine > 50.0 * coarse


def test_abs_is_differentiable_away_from_the_corner():
    assert abs(central_difference(D.absolute, 2.0, D.CORNER_WIDTH) - 1.0) < D.CENTRAL_TOL
    assert abs(central_difference(D.absolute, -2.0, D.CORNER_WIDTH) + 1.0) < D.CENTRAL_TOL


def test_relu_forward_at_zero_is_one():
    assert forward_difference(D.relu, 0.0, D.CORNER_WIDTH) == D.RELU_FORWARD_AT_ZERO


def test_relu_backward_at_zero_is_zero():
    assert backward_difference(D.relu, 0.0, D.CORNER_WIDTH) == D.RELU_BACKWARD_AT_ZERO


def test_relu_central_at_zero_is_one_half():
    assert central_difference(D.relu, 0.0, D.CORNER_WIDTH) == D.RELU_CENTRAL_AT_ZERO


def test_the_relu_central_value_is_neither_defensible_choice():
    value = central_difference(D.relu, 0.0, D.CORNER_WIDTH)
    assert value != 0.0
    assert value != 1.0


@pytest.mark.parametrize(("x", "expected"), [(-1.0, 0.0), (-0.5, 0.0), (0.5, 1.0), (1.0, 1.0)])
def test_relu_is_differentiable_away_from_zero(x, expected):
    assert abs(central_difference(D.relu, x, D.CORNER_WIDTH) - expected) < D.CENTRAL_TOL


def test_relu_on_the_flat_arm_is_exactly_zero():
    assert central_difference(D.relu, -1.0, D.CORNER_WIDTH) == 0.0


# ---------------------------------------------------------------------------
# NumPy
# ---------------------------------------------------------------------------


def test_numpy_gradient_with_scalar_spacing_is_our_central_difference():
    assert numpy_gradient_slope(D.exponential, 1.0, D.COMPARE_WIDTH) == central_difference(
        D.exponential, 1.0, D.COMPARE_WIDTH
    )


@pytest.mark.parametrize("x", [0.5, 1.0, 2.0])
def test_numpy_gradient_matches_at_several_points(x):
    assert numpy_gradient_slope(D.exponential, x, D.COMPARE_WIDTH) == central_difference(
        D.exponential, x, D.COMPARE_WIDTH
    )


def test_numpy_gradient_with_coordinates_differs_in_the_last_bits():
    scalar = numpy_gradient_slope(D.exponential, 1.0, D.COMPARE_WIDTH)
    coords = numpy_gradient_slope_from_coordinates(D.exponential, 1.0, D.COMPARE_WIDTH)
    assert scalar != coords
    assert abs(scalar - coords) < 1e-10


def test_both_numpy_routes_are_within_tolerance_of_the_truth():
    for value in (
        numpy_gradient_slope(D.exponential, 1.0, D.COMPARE_WIDTH),
        numpy_gradient_slope_from_coordinates(D.exponential, 1.0, D.COMPARE_WIDTH),
    ):
        assert abs(value - math.e) < D.CENTRAL_TOL


def test_numpy_gradient_edges_use_a_one_sided_rule():
    h = D.COMPARE_WIDTH
    x = 1.0
    ys = np.array([D.exponential(x - h), D.exponential(x), D.exponential(x + h)])
    edge = float(np.gradient(ys, h)[0])
    assert edge == forward_difference(D.exponential, x - h, h)


def test_numpy_gradient_edges_are_as_bad_as_the_forward_rule():
    h = D.COMPARE_WIDTH
    x = 1.0
    ys = np.array([D.exponential(x - h), D.exponential(x), D.exponential(x + h)])
    grad = np.gradient(ys, h)
    edge_error = abs(float(grad[0]) - math.exp(x - h))
    interior_error = abs(float(grad[1]) - math.e)
    assert edge_error > 1000.0 * interior_error


def test_numpy_is_version_two_or_later():
    assert int(np.__version__.split(".")[0]) >= 2


# ---------------------------------------------------------------------------
# The lab's own honesty
# ---------------------------------------------------------------------------


def test_every_documented_tolerance_is_larger_than_the_error_it_covers():
    assert abs(central_difference(D.exponential, 1.0, D.COMPARE_WIDTH) - math.e) < D.CENTRAL_TOL
    assert abs(forward_difference(D.exponential, 1.0, D.COMPARE_WIDTH) - math.e) < D.FORWARD_TOL
    assert abs(second_difference(D.parabola, 2.0, D.STATIONARY_WIDTH) - 2.0) < D.SECOND_TOL


def test_no_tolerance_is_so_loose_that_it_would_hide_a_real_error():
    # Each tolerance must reject an error one order of magnitude above the
    # bound it was derived from, or it is not testing anything.
    assert D.CENTRAL_TOL < 1e-8
    assert D.FORWARD_TOL < 1e-3
    assert D.SECOND_TOL < 1e-4
    assert D.STATIONARY_TOL < 1e-5


def test_epsilon_is_the_real_float64_epsilon():
    assert D.EPSILON == float(np.finfo(np.float64).eps)


def test_the_settle_widths_shrink_by_a_factor_of_ten_each_time():
    ratios = [D.SETTLE_WIDTHS[i] / D.SETTLE_WIDTHS[i + 1] for i in range(len(D.SETTLE_WIDTHS) - 1)]
    assert all(abs(r - 10.0) < 1e-9 for r in ratios)

lesson_id: D108
day: 108
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-108-derivatives-rates-of-change
  - 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_average_rate_of_change.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_shrinking_intervals.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_rules_checked_numerically.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_forward_and_central.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_the_u_shaped_error.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_zero_derivative_and_curvature.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 07_where_the_derivative_fails.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 -> 97 checks, 0 failure(s), exit 0; pytest examples -> 178 passed; pytest starter -> 1 passed, 99 skipped on an untouched checkout, and 100 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, and the harness was additionally confirmed to exit 0 with no .venv present at all (PYTEST pointing at an interpreter elsewhere), so the clean-disk checks in section 7 cannot fail on a reader who followed the README. Network is needed once to install numpy and pytest; nothing else in the lab opens a socket, and section 7 greps every source file in examples/ and starter/ to prove it. Section 6 re-runs the harness with one expectation deliberately swapped for the belief that the central difference of max(x, 0) at zero is 1.0 rather than 0.5, 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. Measured on this run and REPORTED rather than asserted: the U-shaped error curve over 27 step sizes from 1e-1 to 1e-14 on e**x at x = 1 bottomed out at h = 1.000e-08 with error 6.602751e-09 for the forward difference and at h = 3.162e-06 with error 3.291500e-11 for the central difference, against balance predictions of sqrt(2*eps) = 2.107e-08 and (3*eps)^(1/3) = 8.733e-06 respectively; the tests assert only that each curve is U-shaped with an interior minimum, that both ends are more than a hundred times worse than the middle, that each optimum falls in a documented band, and that each is within a factor of ten of its prediction. Three exact values are asserted with == rather than a tolerance because they are exact in float64 and not approximations: forward_difference(exp, 1.0, 1e-300) is 0.0, because exp(1 + 1e-300) and exp(1) are the same float64 and their difference is exactly zero; the central difference of abs(x) at 0 is 0.0 at every h tried from 1e-2 to 1e-14, where no derivative exists at all; and the central difference of max(x, 0) at 0 is 0.5, which is neither of the two values a framework could defensibly choose. Four findings were discovered while building the lab rather than planned: np.gradient with a scalar spacing is bit-for-bit identical to the from-scratch central difference while np.gradient with an array of the same evenly spaced coordinates differs by 1.879e-11, because the coordinate form takes NumPy''s general unevenly-spaced route; the second difference of abs(x) at zero is exactly 2/h and therefore diverges as h shrinks, which is a far better corner detector than the first derivative''s calm 0.0; the central-difference error curve is not monotone near its minimum, with h = 1e-6 measuring worse than both h = 3.162e-6 and h = 3.162e-7, which is why the U-shape test asks about the ends against the middle rather than for a monotone descent; and x**4 at 0 is a genuine minimum that the second-derivative test cannot distinguish from x**3 at 0, which the suite asserts as "undecided" in both directions rather than glossing. Every tolerance in the lab is derived in examples/dataset.py from the truncation and rounding terms with the arithmetic written out, and a reference test asserts that none of them is loose enough to be meaningless. No deep-learning framework is installed here and no output from one is reproduced anywhere.'

# 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 | Two jobs only. `numpy.gradient` is the library alternative your from-scratch central difference is compared against, bit for bit. And `numpy.finfo(numpy.float64).eps` is where the machine epsilon in `dataset.py` is checked against rather than trusted. |
| `pytest` | 9.1.1 | MIT | Runs both suites: the 178 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.

## Why a version this specific

Less than usual, and it is worth saying so plainly. Almost every number in this
lab comes from `math` and plain float64 arithmetic, and would be identical with
no third-party package at all. NumPy is pinned because two claims depend on it:

- **`np.gradient` with a scalar spacing is bit-for-bit the central difference
  you wrote**, and with an array of coordinates it is not — it takes its general
  unevenly-spaced route and lands a few units in the last place away. Both facts
  are asserted, so a future NumPy that changed either would be reported rather
  than quietly making this lesson wrong.
- **`np.finfo(np.float64).eps` is the value `dataset.EPSILON` claims it is.**
  Every tolerance in the lab is derived from that number, so it is checked
  rather than copied from memory.

The harness also confirms the interpreter's floats are IEEE-754 doubles with a
53-bit significand, because the entire U-shaped error curve is a consequence of
that width.

## 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 most of this lab, which is unusual and is worth taking advantage of.

On a bare `python3` with only the standard library you can write every one of
the ten functions in `starter/derivatives.py`, run the shrinking-interval
sequence, measure the U-shaped error curve across all 27 step sizes, find the
best `h`, classify every stationary point, and reproduce both corner cases. All
of it needs `math` and nothing else.

What you lose: the two `np.gradient` comparisons, the float64 array form of the
error curve, the epsilon cross-check, and the ability to run `pytest`, which
means you would have to call your functions by hand and read the numbers
yourself instead of getting a score.

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

# Watch the Slope Settle — the brief

Seven exercises. Do them in order; each one uses the one before.

Work from the LAB DIRECTORY (the one above this file) and check yourself as
often as you like:

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

On an untouched checkout that prints `1 passed, 99 skipped`. A skip means "not
attempted". A failure means "attempted and wrong", and prints both your answer
and the real one. When it prints `100 passed`, you are finished.

**Do not read `examples/` until you have tried.** The reference is there for
afterwards, and reading it first turns a lab into a transcription exercise.

---

## Exercise 1 — average rates (`derivatives.py`, functions 1.1 and 1.2)

Write `average_rate` and `shrinking_slopes`.

`average_rate(f, a, b)` is rise over run and nothing more. The only interesting
decision is what to do when `a == b`: raise `ZeroDivisionError` with a message
containing the words "interval with width". Do not return 0.0 and do not return
nan. That question genuinely has no answer, and the whole day exists because it
does not.

`shrinking_slopes` is one line built on `average_rate`. Feed it widths that get
smaller and the returned numbers settle. That settling is the limit.

## Exercise 2 — the difference quotients (`derivatives.py`, 2.1 to 2.4)

Write `forward_difference`, `backward_difference`, `central_difference` and
`second_difference`.

Three warnings, all of which have a test named after them:

- The central difference divides by `2 * h`, not by `h`. Forgetting the 2
  doubles every answer you will ever get from it, and doubling is not obviously
  wrong when you do not already know the right value.
- The second difference divides by `h * h`, not by `h`.
- Do not differentiate anything symbolically. These functions may call `f` at
  points you choose and do arithmetic on what comes back, and nothing else.
  That is the honest situation you are in whenever the function is a model
  rather than an equation.

## Exercise 3 — the error curve (`derivatives.py`, 3.1 to 3.3)

Write `error_curve`, `best_step` and `is_u_shaped`.

`error_curve` takes a `rule` — one of your own functions above, passed in as a
value — and returns the absolute error against a slope you already know. It is
only usable when you know the right answer, which is exactly why the lab
measures it on `e**x`: you cannot see the shape of an error you cannot compute.

`is_u_shaped` must be deliberately tolerant. Rounding error near the bottom is a
random walk, not a smooth curve, and a test demanding a monotone descent would
be a test demanding something untrue.

Then run it on the real 27-point grid from `dataset.U_WIDTHS`, which spans
h = 1e-1 down to h = 1e-14, and look at where the bottom is. It is not at the
small end, and understanding why is the point of the day.

## Exercise 4 — flat points (`derivatives.py`, 4.1)

Write `classify_stationary_point`.

Four possible answers, and the fourth is the one that matters: `"undecided"`.
A zero first derivative says the ground is level. It does not say whether you
are at the bottom of a valley, the top of a hill, or on a flat step partway down
a slope. The second derivative resolves two of those three. When it is zero as
well, nothing here can separate `x**3` at 0 — a step — from `x**4` at 0, which
is a genuine minimum.

A test asserts that you return `"undecided"` for `x**4` at 0. Returning
`"minimum"` there would be right by accident, and the suite treats being right
by accident as being wrong.

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

Forty-two predictions, in `answers.py`. Work each one out **before** running
anything. Nearly all can be done on paper.

Two notes on format. Where an option is written across two lines in a comment,
it is one string with single spaces — type it on one line. And where a question
asks for an exception, give the class itself (`ZeroDivisionError`), not a string
naming it.

---

## Order that works

1. Write `average_rate` and `shrinking_slopes`, run the suite, watch nine skips
   turn into passes.
2. Answer exercise 2's predictions while the arithmetic is fresh.
3. Write the three first-difference rules and `second_difference`.
4. Answer exercises 3, 4 and 5.
5. Write `error_curve`, `best_step` and `is_u_shaped`, then run them on the real
   27-point grid and look at the numbers before answering exercise 6.
6. Write `classify_stationary_point` and answer exercise 7.
7. Only now, read `examples/`, run the seven scripts, and see what you would
   have written differently.

## When you are done

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

97 checks, and section 6 of it deliberately breaks one expectation to prove the
harness can go red. A green suite proves nothing until you have watched it fail.

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

Almost all of these can be done on paper or in your head. 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
"""

# =============================================================================
# Exercise 2 -- average rates, before any calculus
# =============================================================================

# 2.1 A car's distance from a post is 4 * t**2 metres at t seconds. What is its
#     average speed, in metres per second, over the whole interval from t = 0
#     to t = 6? A float.
AVERAGE_SPEED_WHOLE_TRIP = None

# 2.2 What is its average speed over the fourth second alone, from t = 3 to
#     t = 4? A float.
AVERAGE_SPEED_FOURTH_SECOND = None

# 2.3 The car was never actually travelling at exactly the answer to 2.1 for
#     the whole trip. What does that average speed describe?
#     One of these strings:
#       "the speed shown on the speedometer at t = 3"
#       "the constant speed that would have covered the same distance in the
#        same time"
#       "the highest speed the car reached"
AVERAGE_SPEED_MEANING = None

# 2.4 What does `average_rate(f, 3.0, 3.0)` do in this lab? Answer with the
#     EXCEPTION CLASS itself, not a string. For example: ValueError
ZERO_WIDTH_RAISES = None

# 2.5 For f(x) = x**2, the average rate over [3, 3 + h] simplifies to a very
#     short expression in h. Which one?
#     One of these strings:  "6", "6 + h", "6 + h**2", "9 + 6h"
SIMPLIFIED_SECANT_SLOPE = None


# =============================================================================
# Exercise 3 -- the shrinking interval and the limit
# =============================================================================

# 3.1 Using your answer to 2.5, what are the four secant slopes for
#     h = 1, 0.1, 0.01 and 0.001? A list of four floats.
SETTLING_SEQUENCE = None

# 3.2 What number does that sequence settle on? A float.
SETTLED_VALUE = None

# 3.3 For f(x) = x**2 at x = 3, do the secant slopes approach 6 from ABOVE or
#     from BELOW as h shrinks through positive values?
#     One of these strings: "above", "below"
APPROACH_DIRECTION = None

# 3.4 The tangent line to y = x**2 at x = 3 passes through (3, 9) with the
#     derivative as its slope. Give (slope, intercept) as a tuple of two
#     floats, for the line written as y = slope * x + intercept.
TANGENT_LINE = None

# 3.5 Which statement about a tangent line is correct?
#     One of these strings:
#       "a tangent line touches the curve at exactly one point"
#       "a tangent line is the line the secant lines approach as the interval
#        shrinks"
#       "a tangent line never crosses the curve"
TANGENT_DEFINITION = None


# =============================================================================
# Exercise 4 -- the rules
# =============================================================================

# 4.1 d/dx of 7, at any x. A float.
DERIVATIVE_OF_SEVEN = None

# 4.2 d/dx of x**5, evaluated at x = 1.5. A float. (Power rule: n * x**(n-1).)
DERIVATIVE_OF_X5_AT_1_5 = None

# 4.3 d/dx of 5 * x**2, evaluated at x = 3. A float.
DERIVATIVE_OF_5X2_AT_3 = None

# 4.4 d/dx of x**2 + x**3, evaluated at x = 2. A float.
DERIVATIVE_OF_SUM_AT_2 = None

# 4.5 d/dx of ln(x), evaluated at x = 4. A float.
DERIVATIVE_OF_LN_AT_4 = None

# 4.6 The slope of b**x at x = 0, for b = 2, is not 1. What number is it?
#     One of these strings:
#       "1, the same for every base"
#       "the natural logarithm of 2, about 0.693"
#       "2, the base itself"
SLOPE_OF_2X_AT_ZERO = None

# 4.7 What makes e special among all the possible bases?
#     One of these strings:
#       "e**x is the only function whose graph is a straight line"
#       "e is the base for which the slope at x = 0 is exactly 1, so e**x is
#        its own derivative"
#       "e is the largest base for which the derivative exists"
WHY_E_IS_SPECIAL = None


# =============================================================================
# Exercise 5 -- forward against central
# =============================================================================

# 5.1 For f(x) = x**2 at x = 3, the forward difference at step h is exactly
#     6 + h. What is the backward difference at the same h?
#     One of these strings: "6 + h", "6 - h", "6", "6 + h**2"
BACKWARD_ON_A_PARABOLA = None

# 5.2 What is the central difference for that same function at that same point,
#     for ANY h at all? A float.
CENTRAL_ON_A_PARABOLA = None

# 5.3 If you divide h by 10, roughly what happens to the FORWARD difference's
#     truncation error?
#     One of these strings: "divided by 10", "divided by 100", "unchanged"
FORWARD_ERROR_SCALING = None

# 5.4 And to the CENTRAL difference's truncation error?
#     One of these strings: "divided by 10", "divided by 100", "unchanged"
CENTRAL_ERROR_SCALING = None

# 5.5 How many calls to f does the central difference need, per estimate?
#     An integer.
CENTRAL_FUNCTION_CALLS = None

# 5.6 The central difference is the average of which two rules?
#     One of these strings:
#       "the forward and backward differences"
#       "the forward difference at h and at 2h"
#       "the first and second differences"
CENTRAL_IS_THE_AVERAGE_OF = None


# =============================================================================
# Exercise 6 -- the U-shaped error curve
# =============================================================================

# 6.1 As h shrinks from 1e-1 towards 1e-14, what does the TRUNCATION error do?
#     One of these strings: "shrinks", "grows", "stays the same"
TRUNCATION_AS_H_SHRINKS = None

# 6.2 And what does the ROUNDING error do?
#     One of these strings: "shrinks", "grows", "stays the same"
#     Note the answer is not the same as 6.1; if it were, there would be no U.
ROUNDING_AS_H_SHRINKS = None

# 6.3 What is `forward_difference(math.exp, 1.0, 1e-300)`? A float, and it is
#     not close to e. Think about what exp(1 + 1e-300) is stored as.
ABSURDLY_SMALL_H_RESULT = None

# 6.4 Why is that the answer?
#     One of these strings:
#       "Python cannot represent 1e-300"
#       "exp(1 + 1e-300) and exp(1) are the same float64, so their difference
#        is exactly zero"
#       "the exponential function is flat near x = 1"
ABSURDLY_SMALL_H_REASON = None

# 6.5 Roughly where does the CENTRAL difference's error bottom out for float64?
#     One of these strings: "around 1e-2", "around 1e-6", "around 1e-16"
BEST_CENTRAL_H_BAND = None

# 6.6 And the FORWARD difference's?
#     One of these strings: "around 1e-2", "around 1e-8", "around 1e-16"
BEST_FORWARD_H_BAND = None

# 6.7 True or false: choosing h = 1e-12 for a central difference is a more
#     careful choice than h = 1e-5. A bool.
TINY_H_IS_MORE_CAREFUL = None


# =============================================================================
# Exercise 7 -- flat points, curvature, and corners
# =============================================================================

# 7.1 f(x) = x**3 - 3x has a zero derivative at x = -1 and at x = +1. Which is
#     the MAXIMUM?
#     One of these strings: "x = -1", "x = +1", "both", "neither"
WHICH_IS_THE_MAXIMUM = None

# 7.2 What is f''(x) for that cubic at x = +1? A float. (f'' of x**3 - 3x is
#     6x, so this is one multiplication.)
SECOND_DERIVATIVE_AT_PLUS_ONE = None

# 7.3 f(x) = x**3 at x = 0 has f'(0) = 0 and f''(0) = 0. What kind of point is
#     it?
#     One of these strings: "minimum", "maximum", "neither"
CUBE_AT_ZERO = None

# 7.4 What does `classify_stationary_point` return there?
#     One of these strings: "minimum", "maximum", "undecided", "not stationary"
CUBE_AT_ZERO_CLASSIFICATION = None

# 7.5 What does a zero first derivative tell you, on its own?
#     One of these strings:
#       "that you are at a minimum"
#       "that the function is flat there, and nothing more"
#       "that the function is constant"
WHAT_ZERO_DERIVATIVE_MEANS = None

# 7.6 f(x) = |x| at x = 0. What does the FORWARD difference return? A float.
ABS_FORWARD_AT_ZERO = None

# 7.7 What does the BACKWARD difference return there? A float.
ABS_BACKWARD_AT_ZERO = None

# 7.8 What does the CENTRAL difference return there? A float.
ABS_CENTRAL_AT_ZERO = None

# 7.9 Does |x| have a derivative at 0? A bool.
ABS_IS_DIFFERENTIABLE_AT_ZERO = None

# 7.10 What does the central difference of max(x, 0) return at x = 0? A float.
RELU_CENTRAL_AT_ZERO = None

# 7.11 The cheapest way to detect that you are standing on a corner, using
#      values you have already computed:
#      One of these strings:
#        "check whether the central difference is zero"
#        "check whether the forward and backward differences disagree"
#        "check whether the function returns nan"
HOW_TO_DETECT_A_CORNER = None

# 7.12 Why does any of this matter for training a model?
#      One of these strings:
#        "the derivative tells you which way to move to make the loss smaller"
#        "the derivative tells you what the loss will be"
#        "the derivative tells you how many layers the network needs"
WHY_DERIVATIVES_MATTER_FOR_AI = None

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

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

"""The invented data, the functions, the step sizes and the tolerances.

Everything in this lab is computed from the definitions below. Nothing is read
from disk, nothing is downloaded, and no number here was chosen to make a test
pass -- every tolerance is derived in `TOLERANCES` from the two error terms that
actually govern a difference quotient, and the derivation is written out beside
it so you can check the arithmetic yourself.

Read this file. Do not change it: the reference tests compare captured values
against the constants here, so editing one moves the goalposts rather than
fixing anything.
"""

from __future__ import annotations

import math

# ---------------------------------------------------------------------------
# The car, which is where the day starts
# ---------------------------------------------------------------------------

# Invented. A car's distance from a marker post, in metres, sampled once a
# second for six seconds. The numbers were chosen so the arithmetic is doable in
# your head: they are 4 * t**2, so the car is accelerating steadily.
CAR_TIMES_S = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
CAR_DISTANCE_M = [0.0, 4.0, 16.0, 36.0, 64.0, 100.0, 144.0]

# Average speed over the whole six seconds: (144 - 0) / (6 - 0).
CAR_AVERAGE_SPEED_WHOLE_TRIP = 24.0
# Average speed over the fourth second only: (64 - 36) / (4 - 3).
CAR_AVERAGE_SPEED_SECOND_FOUR = 28.0
# The speedometer reading at t = 3 exactly, which is d/dt of 4t**2 = 8t.
CAR_INSTANT_SPEED_AT_3 = 24.0


# ---------------------------------------------------------------------------
# The functions the lab differentiates
# ---------------------------------------------------------------------------


def square(x: float) -> float:
    """f(x) = x**2. Exact derivative 2x. The whole lab's first example."""
    return x * x


def square_derivative(x: float) -> float:
    return 2.0 * x


def cubic(x: float) -> float:
    """f(x) = x**3 - 3x. Stationary at x = -1 (a maximum) and x = +1 (a minimum)."""
    return x * x * x - 3.0 * x


def cubic_derivative(x: float) -> float:
    return 3.0 * x * x - 3.0


def cubic_second_derivative(x: float) -> float:
    return 6.0 * x


def parabola(x: float) -> float:
    """f(x) = (x - 2)**2 + 1. Vertex at x = 2, where the slope is exactly zero."""
    return (x - 2.0) ** 2 + 1.0


def parabola_derivative(x: float) -> float:
    return 2.0 * (x - 2.0)


def plain_cube(x: float) -> float:
    """f(x) = x**3. Flat at x = 0 and yet neither a maximum nor a minimum."""
    return x * x * x


def plain_cube_derivative(x: float) -> float:
    return 3.0 * x * x


def exponential(x: float) -> float:
    """f(x) = e**x. Its own derivative, which is what makes e special."""
    return math.exp(x)


def natural_log(x: float) -> float:
    """f(x) = ln(x). Derivative 1/x."""
    return math.log(x)


def absolute(x: float) -> float:
    """f(x) = |x|. No derivative at 0 -- and the corner ReLU is built from."""
    return abs(x)


def relu(x: float) -> float:
    """max(x, 0). Day 102 met this as a transformation; here it is the corner."""
    return x if x > 0.0 else 0.0


# ---------------------------------------------------------------------------
# The points, the widths and the exact answers
# ---------------------------------------------------------------------------

# Where the shrinking-interval demonstration happens.
SETTLE_POINT = 3.0
SETTLE_EXACT_SLOPE = 6.0  # d/dx of x**2 at x = 3
SETTLE_WIDTHS = [1.0, 0.1, 0.01, 0.001]
# For f(x) = x**2 the average rate over [a, a + h] is exactly 2a + h, so at
# a = 3 the sequence below is 6 + h and can be written down without a computer.
SETTLE_EXPECTED_SLOPES = [7.0, 6.1, 6.01, 6.001]

# Where the U-shaped error curve is measured. e**x at x = 1 is chosen because
# every one of its derivatives is e, which makes both error terms easy to state.
U_POINT = 1.0
U_EXACT_SLOPE = math.e  # 2.718281828459045

# 27 step sizes, one per decade and two thirds, spanning 1e-1 down to 1e-14.
U_WIDTHS = [10.0 ** (-1.0 - 0.5 * k) for k in range(27)]

# The point at which forward and central are compared head to head.
COMPARE_WIDTH = 1e-5

# The step used for the stationary-point work. Small enough that the truncation
# term is invisible, large enough that the second difference's h**2 divisor has
# not started amplifying rounding error.
STATIONARY_WIDTH = 1e-4

# Machine epsilon for float64: the gap between 1.0 and the next float up.
EPSILON = 2.220446049250313e-16


# ---------------------------------------------------------------------------
# The tolerances, and where each one comes from
# ---------------------------------------------------------------------------

# A difference quotient carries two errors that pull in opposite directions.
#
#   TRUNCATION comes from the mathematics: the formula is only the limit's
#   approximation at a finite h. Taylor gives, for f = e**x at x = 1 where every
#   derivative equals e:
#       forward:  |error| ~ (h / 2)  * e
#       central:  |error| ~ (h**2 / 6) * e
#
#   ROUNDING comes from the arithmetic: f(x + h) and f(x - h) are each stored to
#   about EPSILON relative precision, their difference cancels most of their
#   digits, and dividing by h magnifies what is left:
#       either rule: |error| ~ EPSILON * e / h
#
# Add the two, put the numbers in, and every tolerance below follows. None was
# reached by running a test and enlarging the number until it went green.

# Central difference of e**x at x = 1 with h = 1e-5:
#   truncation ~ (1e-10 / 6) * e   = 4.5e-11
#   rounding   ~ 2.22e-16 * e / 1e-5 = 6.0e-11
#   sum        ~ 1.1e-10            -> allow 1e-9, roughly nine times the bound
CENTRAL_TOL = 1e-9

# Forward difference of e**x at x = 1 with h = 1e-5:
#   truncation ~ (1e-5 / 2) * e = 1.36e-5
#   rounding   ~ 6.0e-11, negligible beside it
#   sum        ~ 1.4e-5         -> allow 1e-4, seven times the bound
FORWARD_TOL = 1e-4

# The average-rate sequence over [3, 3 + h] for f = x**2. Each term is 2a + h
# computed in float64 from numbers of order 10, so only a few units in the last
# place of about 1e-15 are available to go wrong.
EXACT_TOL = 1e-12

# The second difference of a cubic is exact in the mathematics -- the h**2 term
# in its Taylor expansion is the answer and there is no h**4 term to leave
# behind -- so only rounding is in play. With h = 1e-4 and |f| of order 2:
#   rounding ~ 4 * EPSILON * |f| / h**2 = 4 * 2.22e-16 * 2 / 1e-8 = 1.8e-7
#   -> allow 1e-5, about fifty times the bound
SECOND_TOL = 1e-5

# A stationary point found with the central difference at h = 1e-4 on the cubic:
# the truncation term is exactly h**2 * f'''/6 = h**2, which is 1e-8.
#   -> allow 1e-6, a hundred times the bound
STATIONARY_TOL = 1e-6


# ---------------------------------------------------------------------------
# The derivative rules, stated as facts and checked numerically in the lab
# ---------------------------------------------------------------------------

# (name, f, exact f', a point to check it at)
RULE_CASES = [
    ("constant: d/dx of 7 is 0", lambda x: 7.0, lambda x: 0.0, 2.0),
    ("power: d/dx of x**2 is 2x", square, square_derivative, 3.0),
    ("power: d/dx of x**5 is 5x**4", lambda x: x**5, lambda x: 5.0 * x**4, 1.5),
    ("power: d/dx of 1/x is -1/x**2", lambda x: 1.0 / x, lambda x: -1.0 / (x * x), 2.0),
    ("constant multiple: d/dx of 5x**2 is 10x", lambda x: 5.0 * x * x, lambda x: 10.0 * x, 3.0),
    ("sum: d/dx of x**2 + x**3 is 2x + 3x**2", lambda x: x**2 + x**3, lambda x: 2.0 * x + 3.0 * x**2, 2.0),
    ("exponential: d/dx of e**x is e**x", exponential, exponential, 1.0),
    ("logarithm: d/dx of ln(x) is 1/x", natural_log, lambda x: 1.0 / x, 4.0),
]

# The exact slopes those eight cases must produce, written out so a reader can
# check them by hand rather than by rerunning the lab.
RULE_EXPECTED = [0.0, 6.0, 25.3125, -0.25, 30.0, 16.0, math.e, 0.25]

# The eight rule cases are checked with the central difference at h = 1e-5. The
# worst truncation term among them belongs to x**5 at 1.5, where the third
# derivative is 60 * 1.5**2 = 135:
#   truncation ~ (h**2 / 6) * 135 = (1e-10 / 6) * 135 = 2.25e-9
#   -> allow 1e-8, about four times the bound
RULE_TOL = 1e-8

# The corner cases at x = 0, where no derivative exists.
CORNER_WIDTH = 1e-5
# |x|: the one-sided slopes are -1 and +1 and disagree, which is the whole
# reason there is no derivative. The central difference averages them to zero
# and reports that zero with total confidence.
ABS_FORWARD_AT_ZERO = 1.0
ABS_BACKWARD_AT_ZERO = -1.0
ABS_CENTRAL_AT_ZERO = 0.0
# max(x, 0): the one-sided slopes are 0 and 1, so the central difference gives
# their average, 0.5. Deep-learning frameworks do not use 0.5; they pick one of
# the one-sided values and move on, and that choice is a convention rather than
# a theorem.
RELU_FORWARD_AT_ZERO = 1.0
RELU_BACKWARD_AT_ZERO = 0.0
RELU_CENTRAL_AT_ZERO = 0.5

"""Exercise 1 -- ten functions to write. Your work goes here.

Each function raises NotImplementedError until you write it, and the test suite
SKIPS anything still unwritten rather than failing it. So 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 you write the body. Every one gives you the formula
in words, and a worked example small enough to check on paper.

One rule for the whole file: no calculus. Nothing here differentiates a formula
symbolically. Every function is allowed to do exactly one thing -- call `f` at
points you choose and do arithmetic on what comes back -- which is the honest
situation you are in whenever the function is a model rather than an equation.
"""

from __future__ import annotations

from collections.abc import Callable, Sequence

import numpy as np

Function = Callable[[float], float]


# ===========================================================================
# 1. Average rate of change
# ===========================================================================


def average_rate(f: Function, a: float, b: float) -> float:
    """1.1 -- rise over run between a and b.

    How much the output changed, divided by how much the input changed:

        ( f(b) - f(a) ) / ( b - a )

    If a and b are equal there is no interval and no answer, so RAISE
    ZeroDivisionError with a message containing the words "interval with
    width". Do not return 0.0 and do not return nan: the point of the whole day
    is that this question genuinely has no answer, and a function that invents
    one is teaching the wrong lesson.

    >>> average_rate(lambda x: x * x, 3.0, 4.0)
    7.0
    >>> average_rate(lambda x: 4.0 * x * x, 0.0, 6.0)
    24.0
    """
    raise NotImplementedError("average_rate")


def shrinking_slopes(f: Function, a: float, widths: Sequence[float]) -> list[float]:
    """1.2 -- the average rate over [a, a + h], one entry per h, in order.

    One line if you use a list comprehension over `widths` and call
    `average_rate`. Do not reimplement rise over run here.

    >>> shrinking_slopes(lambda x: x * x, 3.0, [1.0, 0.1])
    [7.0, 6.100000000000007]
    """
    raise NotImplementedError("shrinking_slopes")


# ===========================================================================
# 2. The three difference quotients
# ===========================================================================


def forward_difference(f: Function, x: float, h: float) -> float:
    """2.1 -- the slope of the secant from x to x + h.

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

    >>> forward_difference(lambda x: x * x, 3.0, 0.1)
    6.100000000000012
    """
    raise NotImplementedError("forward_difference")


def backward_difference(f: Function, x: float, h: float) -> float:
    """2.2 -- the slope of the secant from x - h to x.

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

    >>> backward_difference(lambda x: x * x, 3.0, 0.1)
    5.899999999999999
    """
    raise NotImplementedError("backward_difference")


def central_difference(f: Function, x: float, h: float) -> float:
    """2.3 -- the slope of the secant from x - h to x + h, straddling x.

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

    Watch the 2. Forgetting it halves every answer you get, and halving is not
    obviously wrong when you do not already know the right value -- which is
    the single most common bug in this whole topic.

    >>> central_difference(lambda x: x * x, 3.0, 0.1)
    6.000000000000005
    """
    raise NotImplementedError("central_difference")


def second_difference(f: Function, x: float, h: float) -> float:
    """2.4 -- an estimate of f''(x): the rate of change of the rate of change.

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

    Read it as: how far does the middle value sag below the average of its two
    neighbours? Sagging down is positive, a bowl. Bulging up is negative, a
    dome.

    Note the h SQUARED in the divisor, which is why this rule is far more
    sensitive to a badly chosen h than the first-difference rules are.

    >>> round(second_difference(lambda x: x * x, 3.0, 0.001), 6)
    2.0
    """
    raise NotImplementedError("second_difference")


# ===========================================================================
# 3. Measuring the error
# ===========================================================================


def error_curve(
    f: Function,
    x: float,
    exact_slope: float,
    widths: Sequence[float],
    rule: Callable[[Function, float, float], float],
) -> list[float]:
    """3.1 -- the absolute error of `rule` against a known answer, one per width.

    `rule` is one of your own functions above, passed in as a value. Call it as
    `rule(f, x, h)`.

    >>> error_curve(lambda x: x * x, 3.0, 6.0, [1.0, 0.1], forward_difference)
    [1.0, 0.10000000000001208]
    """
    raise NotImplementedError("error_curve")


def best_step(widths: Sequence[float], errors: Sequence[float]) -> tuple[float, float]:
    """3.2 -- the (width, error) pair with the smallest error: the bottom of the U.

    Raise ValueError if the two sequences differ in length, and ValueError if
    they are empty. Ties go to the FIRST minimum found, which for a descending
    list of widths is the largest h that achieves it -- the safer choice,
    because a larger step sits further from the cancellation cliff.

    >>> best_step([1.0, 0.1, 0.01], [5.0, 0.5, 2.0])
    (0.1, 0.5)
    """
    raise NotImplementedError("best_step")


def is_u_shaped(errors: Sequence[float]) -> bool:
    """3.3 -- True when the errors fall to an interior minimum and rise again.

    Be deliberately tolerant. Rounding error is a random walk rather than a
    smooth curve, so do NOT demand a monotone descent -- that would be a test
    demanding something untrue, which is worse than no test.

    Ask three things:
      * there are at least three entries;
      * the smallest error is NOT the first or the last entry;
      * both the first and the last entry are more than ten times the smallest.

    >>> is_u_shaped([100.0, 1.0, 0.01, 1.0, 100.0])
    True
    >>> is_u_shaped([100.0, 10.0, 1.0, 0.1])
    False
    """
    raise NotImplementedError("is_u_shaped")


# ===========================================================================
# 4. What a zero derivative does and does not tell you
# ===========================================================================


def classify_stationary_point(f: Function, x: float, h: float, tol: float) -> str:
    """4.1 -- name the kind of point x is, from the first two derivatives.

    Return exactly one of these strings:

        "not stationary"  the central difference is larger than tol in size
        "minimum"         flat, and the second difference is above  +tol
        "maximum"         flat, and the second difference is below  -tol
        "undecided"       flat, and the second difference is inside +/- tol

    That last one is the important one, and returning it is not a cop-out. A
    zero first derivative says the ground is level. It does not say whether you
    are at the bottom of a valley, the top of a hill, or on a flat step partway
    down a slope, and when the curvature is zero as well, nothing here can
    separate x**3 at 0 (a step) from x**4 at 0 (a genuine minimum). Reporting
    "minimum" there would be a confident lie.

    >>> classify_stationary_point(lambda x: (x - 2.0) ** 2 + 1.0, 2.0, 1e-4, 1e-6)
    'minimum'
    >>> classify_stationary_point(lambda x: x ** 3, 0.0, 1e-4, 1e-6)
    'undecided'
    """
    raise NotImplementedError("classify_stationary_point")


# ===========================================================================
# Written for you -- read these, the tests use them
# ===========================================================================


def tangent_at(f: Function, x: float, h: float) -> tuple[float, float]:
    """(slope, intercept) of the tangent line to f at x.

    The tangent passes through (x, f(x)) with the derivative as its slope, so
    rearranging y = mx + c gives c = f(x) - m*x. This calls YOUR
    central_difference, so it starts working the moment exercise 2.3 does.
    """
    slope = central_difference(f, x, h)
    return (slope, f(x) - slope * x)


def numpy_error_curve(
    f: Function,
    x: float,
    exact_slope: float,
    widths: Sequence[float],
    rule: Callable[[Function, float, float], float],
) -> np.ndarray:
    """Your error_curve as a float64 array, so `.argmin()` and a plot are one line."""
    return np.array(error_curve(f, x, exact_slope, widths, rule), dtype=np.float64)

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

    .venv/bin/pytest starter -q

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

Every test that exercises your code runs its whole body inside `written(...)`,
so a test skips 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. That would say
"attempted and wrong" about work you had not attempted, which is precisely the
lie this suite exists to avoid.
"""

from __future__ import annotations

import math

import numpy as np
import pytest

import answers
import dataset
from derivatives import (
    average_rate,
    backward_difference,
    best_step,
    central_difference,
    classify_stationary_point,
    error_curve,
    forward_difference,
    is_u_shaped,
    numpy_error_curve,
    second_difference,
    shrinking_slopes,
    tangent_at,
)


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


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


def car_distance(t: float) -> float:
    return 4.0 * t * t


# -- 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 dataset.SETTLE_EXPECTED_SLOPES == [7.0, 6.1, 6.01, 6.001], "dataset.py loads"
    assert dataset.square(3.0) == 9.0, "the functions load"


# -- Exercise 1.1: average_rate ----------------------------------------------


def test_1_1_average_rate_over_the_whole_trip():
    assert written(average_rate, car_distance, 0.0, 6.0) == 24.0


def test_1_1_average_rate_over_the_fourth_second():
    assert written(average_rate, car_distance, 3.0, 4.0) == 28.0


def test_1_1_average_rate_of_a_straight_line_is_its_slope():
    def line(x: float) -> float:
        return 3.0 * x + 5.0

    assert written(average_rate, line, -100.0, 100.0) == 3.0


def test_1_1_average_rate_does_not_care_which_end_you_start_at():
    forwards = written(average_rate, dataset.square, 1.0, 4.0)
    backwards = written(average_rate, dataset.square, 4.0, 1.0)
    assert forwards == backwards == 5.0


def test_1_1_average_rate_refuses_a_zero_width_interval():
    def call():
        with pytest.raises(ZeroDivisionError) as caught:
            average_rate(car_distance, 3.0, 3.0)
        return str(caught.value)

    message = written(call)
    assert "interval with width" in message, (
        "raise ZeroDivisionError with a message containing 'interval with width'"
    )


# -- Exercise 1.2: shrinking_slopes ------------------------------------------


def test_1_2_shrinking_slopes_returns_one_value_per_width():
    slopes = written(shrinking_slopes, dataset.square, 3.0, dataset.SETTLE_WIDTHS)
    assert len(slopes) == len(dataset.SETTLE_WIDTHS)


def test_1_2_shrinking_slopes_are_six_plus_h():
    slopes = written(shrinking_slopes, dataset.square, 3.0, dataset.SETTLE_WIDTHS)
    for slope, expected in zip(slopes, dataset.SETTLE_EXPECTED_SLOPES):
        assert abs(slope - expected) < dataset.EXACT_TOL


def test_1_2_shrinking_slopes_settle_on_the_derivative():
    slopes = written(shrinking_slopes, dataset.square, 3.0, dataset.SETTLE_WIDTHS)
    gaps = [abs(s - 6.0) for s in slopes]
    assert gaps == sorted(gaps, reverse=True), "each interval must be closer than the last"
    assert gaps[-1] < 0.002


def test_1_2_shrinking_slopes_keeps_the_order_it_was_given():
    slopes = written(shrinking_slopes, dataset.square, 3.0, [0.001, 1.0])
    assert slopes[0] < slopes[1], "return the results in the order the widths came in"


# -- Exercise 2.1 to 2.3: the three difference quotients ---------------------


@pytest.mark.parametrize("h", [1.0, 0.1, 0.01, 0.001])
def test_2_1_forward_difference_of_a_parabola(h):
    assert abs(written(forward_difference, dataset.square, 3.0, h) - (6.0 + h)) < 1e-9


@pytest.mark.parametrize("h", [1.0, 0.1, 0.01, 0.001])
def test_2_2_backward_difference_of_a_parabola(h):
    assert abs(written(backward_difference, dataset.square, 3.0, h) - (6.0 - h)) < 1e-9


@pytest.mark.parametrize("h", [1.0, 0.1, 0.01, 0.001])
def test_2_3_central_difference_is_exact_on_a_parabola(h):
    assert abs(written(central_difference, dataset.square, 3.0, h) - 6.0) < 1e-11


def test_2_3_central_difference_did_not_forget_the_two():
    """The commonest bug in the topic: dividing by h rather than by 2h."""
    value = written(central_difference, dataset.square, 3.0, 0.1)
    assert abs(value - 12.0) > 1.0, "this is exactly double the right answer: divide by 2 * h"
    assert abs(value - 6.0) < 1e-11


def test_2_3_central_difference_is_the_average_of_the_other_two():
    h = 0.1
    forward = written(forward_difference, dataset.exponential, 1.0, h)
    backward = written(backward_difference, dataset.exponential, 1.0, h)
    central = written(central_difference, dataset.exponential, 1.0, h)
    assert abs((forward + backward) / 2.0 - central) < 1e-14


def test_2_3_central_meets_its_documented_tolerance_on_the_exponential():
    error = abs(written(central_difference, dataset.exponential, 1.0, dataset.COMPARE_WIDTH) - math.e)
    assert error < dataset.CENTRAL_TOL


def test_2_1_forward_meets_its_documented_tolerance_on_the_exponential():
    error = abs(written(forward_difference, dataset.exponential, 1.0, dataset.COMPARE_WIDTH) - math.e)
    assert error < dataset.FORWARD_TOL


def test_2_3_central_beats_forward_by_a_thousandfold_at_the_same_step():
    h = dataset.COMPARE_WIDTH
    forward_error = abs(written(forward_difference, dataset.exponential, 1.0, h) - math.e)
    central_error = abs(written(central_difference, dataset.exponential, 1.0, h) - math.e)
    assert central_error * 1000.0 < forward_error


# -- Exercise 2.4: second_difference -----------------------------------------


def test_2_4_second_difference_of_a_parabola_is_two():
    value = written(second_difference, dataset.parabola, 2.0, dataset.STATIONARY_WIDTH)
    assert abs(value - 2.0) < dataset.SECOND_TOL


def test_2_4_second_difference_at_the_cubic_minimum_is_six():
    value = written(second_difference, dataset.cubic, 1.0, dataset.STATIONARY_WIDTH)
    assert abs(value - 6.0) < dataset.SECOND_TOL


def test_2_4_second_difference_at_the_cubic_maximum_is_minus_six():
    value = written(second_difference, dataset.cubic, -1.0, dataset.STATIONARY_WIDTH)
    assert abs(value + 6.0) < dataset.SECOND_TOL


def test_2_4_second_difference_did_not_forget_to_square_h():
    value = written(second_difference, dataset.parabola, 2.0, 0.01)
    assert abs(value - 0.02) > 1e-3, "dividing by h rather than h**2 gives 0.02 here"
    assert abs(value - 2.0) < dataset.SECOND_TOL


def test_2_4_second_difference_tracks_the_exact_second_derivative():
    for x in [-2.0, -0.5, 0.5, 2.0]:
        value = written(second_difference, dataset.cubic, x, dataset.STATIONARY_WIDTH)
        assert abs(value - dataset.cubic_second_derivative(x)) < dataset.SECOND_TOL


# -- Exercise 3.1: error_curve -----------------------------------------------


def test_3_1_error_curve_returns_one_error_per_width():
    errors = written(error_curve, dataset.exponential, 1.0, math.e, dataset.U_WIDTHS, central_difference)
    assert len(errors) == len(dataset.U_WIDTHS)


def test_3_1_error_curve_values_are_non_negative():
    errors = written(error_curve, dataset.exponential, 1.0, math.e, dataset.U_WIDTHS, central_difference)
    assert all(e >= 0.0 for e in errors), "an absolute error is never negative"


def test_3_1_error_curve_uses_the_rule_it_was_handed():
    forward_errors = written(error_curve, dataset.exponential, 1.0, math.e, [1e-3], forward_difference)
    central_errors = written(error_curve, dataset.exponential, 1.0, math.e, [1e-3], central_difference)
    assert forward_errors[0] > central_errors[0] * 100.0


def test_3_1_error_curve_on_a_parabola_is_exactly_h_for_the_forward_rule():
    errors = written(error_curve, dataset.square, 3.0, 6.0, [1.0, 0.1, 0.01], forward_difference)
    for error, h in zip(errors, [1.0, 0.1, 0.01]):
        assert abs(error - h) < 1e-9


# -- Exercise 3.2: best_step -------------------------------------------------


def test_3_2_best_step_finds_the_smallest_error():
    assert written(best_step, [1.0, 0.1, 0.01], [5.0, 0.5, 2.0]) == (0.1, 0.5)


def test_3_2_best_step_prefers_the_first_of_a_tie():
    assert written(best_step, [1.0, 0.1, 0.01], [2.0, 0.5, 0.5]) == (0.1, 0.5)


def test_3_2_best_step_rejects_mismatched_lengths():
    def call():
        with pytest.raises(ValueError):
            best_step([1.0, 2.0], [1.0])
        return True

    assert written(call)


def test_3_2_best_step_rejects_an_empty_grid():
    def call():
        with pytest.raises(ValueError):
            best_step([], [])
        return True

    assert written(call)


def test_3_2_best_step_on_the_real_central_curve():
    errors = written(error_curve, dataset.exponential, 1.0, math.e, dataset.U_WIDTHS, central_difference)
    h, error = written(best_step, dataset.U_WIDTHS, errors)
    assert error == min(errors)
    assert 1e-7 <= h <= 1e-4, "the bottom of the U for a central difference in float64"


def test_3_2_best_step_on_the_real_forward_curve():
    errors = written(error_curve, dataset.exponential, 1.0, math.e, dataset.U_WIDTHS, forward_difference)
    h, _ = written(best_step, dataset.U_WIDTHS, errors)
    assert 1e-9 <= h <= 1e-6, "the bottom of the U for a forward difference in float64"


# -- Exercise 3.3: is_u_shaped -----------------------------------------------


def test_3_3_is_u_shaped_accepts_a_clear_u():
    assert written(is_u_shaped, [100.0, 1.0, 0.01, 1.0, 100.0]) is True


def test_3_3_is_u_shaped_rejects_a_monotone_fall():
    assert written(is_u_shaped, [100.0, 10.0, 1.0, 0.1]) is False


def test_3_3_is_u_shaped_rejects_a_monotone_rise():
    assert written(is_u_shaped, [0.1, 1.0, 10.0, 100.0]) is False


def test_3_3_is_u_shaped_rejects_a_sequence_that_is_too_short():
    assert written(is_u_shaped, [1.0, 0.5]) is False


def test_3_3_the_real_central_curve_is_u_shaped():
    errors = written(error_curve, dataset.exponential, 1.0, math.e, dataset.U_WIDTHS, central_difference)
    assert written(is_u_shaped, errors) is True


def test_3_3_the_real_forward_curve_is_u_shaped():
    errors = written(error_curve, dataset.exponential, 1.0, math.e, dataset.U_WIDTHS, forward_difference)
    assert written(is_u_shaped, errors) is True


def test_3_3_both_ends_of_the_real_curve_are_far_worse_than_the_middle():
    errors = written(error_curve, dataset.exponential, 1.0, math.e, dataset.U_WIDTHS, central_difference)
    best = min(errors)
    assert errors[0] > 100.0 * best
    assert errors[-1] > 100.0 * best


# -- Exercise 4.1: classify_stationary_point ---------------------------------


@pytest.mark.parametrize(
    ("f", "x", "expected"),
    [
        (dataset.parabola, 2.0, "minimum"),
        (dataset.cubic, 1.0, "minimum"),
        (dataset.cubic, -1.0, "maximum"),
        (dataset.plain_cube, 0.0, "undecided"),
        (dataset.cubic, 0.0, "not stationary"),
        (dataset.exponential, 0.0, "not stationary"),
    ],
    ids=["parabola vertex", "cubic minimum", "cubic maximum", "flat step", "sloping cubic", "exponential"],
)
def test_4_1_classification(f, x, expected):
    verdict = written(classify_stationary_point, f, x, dataset.STATIONARY_WIDTH, dataset.STATIONARY_TOL)
    assert verdict == expected


def test_4_1_a_genuine_minimum_with_zero_curvature_is_also_undecided():
    verdict = written(
        classify_stationary_point, lambda x: x**4, 0.0, dataset.STATIONARY_WIDTH, dataset.STATIONARY_TOL
    )
    assert verdict == "undecided", (
        "x**4 at 0 IS a minimum, and no second derivative can show it -- "
        "reporting 'minimum' here would be a lie you got away with"
    )


# -- The written-for-you helpers, which start working when your code does ----


def test_tangent_at_three_is_six_x_minus_nine():
    slope, intercept = written(tangent_at, dataset.square, 3.0, dataset.COMPARE_WIDTH)
    assert abs(slope - 6.0) < dataset.CENTRAL_TOL
    assert abs(intercept + 9.0) < 1e-8


def test_numpy_error_curve_is_a_float64_array_of_your_errors():
    array = written(
        numpy_error_curve, dataset.exponential, 1.0, math.e, dataset.U_WIDTHS, central_difference
    )
    assert array.dtype == np.float64
    assert array.shape == (len(dataset.U_WIDTHS),)
    assert 1e-7 <= dataset.U_WIDTHS[int(array.argmin())] <= 1e-4


# -- Exercise 2: average rates ------------------------------------------------


def test_2_1_prediction_average_speed_whole_trip():
    assert predicted("AVERAGE_SPEED_WHOLE_TRIP") == 24.0


def test_2_2_prediction_average_speed_fourth_second():
    assert predicted("AVERAGE_SPEED_FOURTH_SECOND") == 28.0


def test_2_3_prediction_what_an_average_speed_means():
    assert predicted("AVERAGE_SPEED_MEANING") == (
        "the constant speed that would have covered the same distance in the same time"
    )


def test_2_4_prediction_zero_width_raises():
    assert predicted("ZERO_WIDTH_RAISES") is ZeroDivisionError


def test_2_5_prediction_simplified_secant_slope():
    assert predicted("SIMPLIFIED_SECANT_SLOPE") == "6 + h"


# -- Exercise 3: the shrinking interval ---------------------------------------


def test_3_1_prediction_settling_sequence():
    assert predicted("SETTLING_SEQUENCE") == [7.0, 6.1, 6.01, 6.001]


def test_3_2_prediction_settled_value():
    assert predicted("SETTLED_VALUE") == 6.0


def test_3_3_prediction_approach_direction():
    assert predicted("APPROACH_DIRECTION") == "above"


def test_3_4_prediction_tangent_line():
    assert predicted("TANGENT_LINE") == (6.0, -9.0)


def test_3_5_prediction_tangent_definition():
    assert predicted("TANGENT_DEFINITION") == (
        "a tangent line is the line the secant lines approach as the interval shrinks"
    )


# -- Exercise 4: the rules -----------------------------------------------------


def test_4_1_prediction_derivative_of_seven():
    assert predicted("DERIVATIVE_OF_SEVEN") == 0.0


def test_4_2_prediction_derivative_of_x5():
    assert predicted("DERIVATIVE_OF_X5_AT_1_5") == 25.3125


def test_4_3_prediction_derivative_of_5x2():
    assert predicted("DERIVATIVE_OF_5X2_AT_3") == 30.0


def test_4_4_prediction_derivative_of_sum():
    assert predicted("DERIVATIVE_OF_SUM_AT_2") == 16.0


def test_4_5_prediction_derivative_of_ln():
    assert predicted("DERIVATIVE_OF_LN_AT_4") == 0.25


def test_4_6_prediction_slope_of_2x_at_zero():
    assert predicted("SLOPE_OF_2X_AT_ZERO") == "the natural logarithm of 2, about 0.693"


def test_4_7_prediction_why_e_is_special():
    assert predicted("WHY_E_IS_SPECIAL") == (
        "e is the base for which the slope at x = 0 is exactly 1, so e**x is its own derivative"
    )


# -- Exercise 5: forward against central ---------------------------------------


def test_5_1_prediction_backward_on_a_parabola():
    assert predicted("BACKWARD_ON_A_PARABOLA") == "6 - h"


def test_5_2_prediction_central_on_a_parabola():
    assert predicted("CENTRAL_ON_A_PARABOLA") == 6.0


def test_5_3_prediction_forward_error_scaling():
    assert predicted("FORWARD_ERROR_SCALING") == "divided by 10"


def test_5_4_prediction_central_error_scaling():
    assert predicted("CENTRAL_ERROR_SCALING") == "divided by 100"


def test_5_5_prediction_central_function_calls():
    assert predicted("CENTRAL_FUNCTION_CALLS") == 2


def test_5_6_prediction_central_is_the_average_of():
    assert predicted("CENTRAL_IS_THE_AVERAGE_OF") == "the forward and backward differences"


# -- Exercise 6: the U ---------------------------------------------------------


def test_6_1_prediction_truncation_as_h_shrinks():
    assert predicted("TRUNCATION_AS_H_SHRINKS") == "shrinks"


def test_6_2_prediction_rounding_as_h_shrinks():
    assert predicted("ROUNDING_AS_H_SHRINKS") == "grows"


def test_6_3_prediction_absurdly_small_h_result():
    assert predicted("ABSURDLY_SMALL_H_RESULT") == 0.0


def test_6_4_prediction_absurdly_small_h_reason():
    assert predicted("ABSURDLY_SMALL_H_REASON") == (
        "exp(1 + 1e-300) and exp(1) are the same float64, so their difference is exactly zero"
    )


def test_6_5_prediction_best_central_h_band():
    assert predicted("BEST_CENTRAL_H_BAND") == "around 1e-6"


def test_6_6_prediction_best_forward_h_band():
    assert predicted("BEST_FORWARD_H_BAND") == "around 1e-8"


def test_6_7_prediction_tiny_h_is_not_more_careful():
    assert predicted("TINY_H_IS_MORE_CAREFUL") is False


# -- Exercise 7: flat points and corners ---------------------------------------


def test_7_1_prediction_which_is_the_maximum():
    assert predicted("WHICH_IS_THE_MAXIMUM") == "x = -1"


def test_7_2_prediction_second_derivative_at_plus_one():
    assert predicted("SECOND_DERIVATIVE_AT_PLUS_ONE") == 6.0


def test_7_3_prediction_cube_at_zero():
    assert predicted("CUBE_AT_ZERO") == "neither"


def test_7_4_prediction_cube_at_zero_classification():
    assert predicted("CUBE_AT_ZERO_CLASSIFICATION") == "undecided"


def test_7_5_prediction_what_zero_derivative_means():
    assert predicted("WHAT_ZERO_DERIVATIVE_MEANS") == "that the function is flat there, and nothing more"


def test_7_6_prediction_abs_forward_at_zero():
    assert predicted("ABS_FORWARD_AT_ZERO") == 1.0


def test_7_7_prediction_abs_backward_at_zero():
    assert predicted("ABS_BACKWARD_AT_ZERO") == -1.0


def test_7_8_prediction_abs_central_at_zero():
    assert predicted("ABS_CENTRAL_AT_ZERO") == 0.0


def test_7_9_prediction_abs_is_not_differentiable_at_zero():
    assert predicted("ABS_IS_DIFFERENTIABLE_AT_ZERO") is False


def test_7_10_prediction_relu_central_at_zero():
    assert predicted("RELU_CENTRAL_AT_ZERO") == 0.5


def test_7_11_prediction_how_to_detect_a_corner():
    assert predicted("HOW_TO_DETECT_A_CORNER") == "check whether the forward and backward differences disagree"


def test_7_12_prediction_why_derivatives_matter_for_ai():
    assert predicted("WHY_DERIVATIVES_MATTER_FOR_AI") == (
        "the derivative tells you which way to move to make the loss smaller"
    )

#!/usr/bin/env bash
# Tests for the Day 108 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:
#
#   * the average rate of x**2 over [3, 3 + h] is 6 + h for every h, and the
#     four-term sequence 7, 6.1, 6.01, 6.001 gets closer to 6 every time;
#   * the central difference is exact on a parabola and beats the forward
#     difference by more than two hundred thousand times on e**x at h = 1e-5;
#   * the error across 27 step sizes from 1e-1 to 1e-14 is U-shaped -- it
#     falls, reaches an interior minimum, and rises again -- and the measured
#     best h is reported rather than asserted to a fixed value;
#   * a forward difference at h = 1e-300 returns exactly 0.0, with no warning;
#   * the first derivative is zero at a minimum, a maximum and a flat step
#     alike, and only the second derivative tells them apart -- and at x**3
#     and x**4 it cannot, which the suite asserts rather than glosses;
#   * the central difference of |x| at 0 returns 0.0 and the one-sided rules
#     return +1 and -1, so a value was produced where no derivative exists;
#   * 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 108 — Watch the Slope Settle"
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}"

float_width="$("${python_bin}" -c "import sys; print(sys.float_info.mant_dig)")"
check_eq "Python floats are IEEE-754 doubles with a 53-bit significand" "53" "${float_width}"

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

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

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

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

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

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

import numpy as np

import dataset as D
from derivatives import (
    average_rate,
    backward_difference,
    best_step,
    central_difference,
    classify_stationary_point,
    error_curve,
    forward_difference,
    is_u_shaped,
    numpy_gradient_slope,
    numpy_gradient_slope_from_coordinates,
    second_difference,
    shrinking_slopes,
    tangent_at,
)


def car(t):
    return 4.0 * t * t


# -- average rates ---------------------------------------------------------
print("car_whole_trip", average_rate(car, 0.0, 6.0))
print("car_fourth_second", average_rate(car, 3.0, 4.0))
print("car_third_second", average_rate(car, 2.0, 3.0))
try:
    average_rate(car, 3.0, 3.0)
except Exception as exc:  # deliberately broad: the TYPE is what is asserted
    print("zero_width_raises", type(exc).__name__)
else:
    print("zero_width_raises", "NOTHING_RAISED")

# -- the shrinking sequence -------------------------------------------------
slopes = shrinking_slopes(D.square, D.SETTLE_POINT, D.SETTLE_WIDTHS)
print("settle_sequence", "|".join(f"{s:.6f}" for s in slopes))
print("settle_matches_six_plus_h",
      all(abs(s - e) < D.EXACT_TOL for s, e in zip(slopes, D.SETTLE_EXPECTED_SLOPES)))
gaps = [abs(s - 6.0) for s in slopes]
print("settle_monotone", gaps == sorted(gaps, reverse=True))
print("settle_from_above", all(s > 6.0 for s in slopes))
left = shrinking_slopes(D.square, D.SETTLE_POINT, [-w for w in D.SETTLE_WIDTHS])
print("settle_from_below", all(s < 6.0 for s in left))
slope, intercept = tangent_at(D.square, 3.0, D.COMPARE_WIDTH)
print("tangent_slope", round(slope, 6))
print("tangent_intercept", round(intercept, 6))

# -- the rules --------------------------------------------------------------
measured = [central_difference(f, x, D.COMPARE_WIDTH) for _, f, _, x in D.RULE_CASES]
print("rules_all_within_tolerance",
      all(abs(m - e) < D.RULE_TOL for m, e in zip(measured, D.RULE_EXPECTED)))
print("rules_exact_values", "|".join(f"{v:.4f}" for v in D.RULE_EXPECTED))
print("slope_of_two_to_the_x_at_zero", round(central_difference(lambda x: 2.0**x, 0.0, D.COMPARE_WIDTH), 9))
print("ln_two", round(math.log(2.0), 9))
print("slope_of_e_to_the_x_at_zero", round(central_difference(D.exponential, 0.0, D.COMPARE_WIDTH), 9))

# -- forward, backward, central ---------------------------------------------
print("forward_on_parabola", forward_difference(D.square, 3.0, 0.1))
print("backward_on_parabola", backward_difference(D.square, 3.0, 0.1))
print("central_on_parabola_is_six", abs(central_difference(D.square, 3.0, 0.1) - 6.0) < 1e-11)
fwd_err = abs(forward_difference(D.exponential, 1.0, D.COMPARE_WIDTH) - math.e)
cen_err = abs(central_difference(D.exponential, 1.0, D.COMPARE_WIDTH) - math.e)
print("forward_error_at_1e5", f"{fwd_err:.6e}")
print("central_error_at_1e5", f"{cen_err:.6e}")
print("central_beats_forward_thousandfold", cen_err * 1000.0 < fwd_err)
print("central_advantage_ratio", int(fwd_err / cen_err))
coarse = abs(central_difference(D.exponential, 1.0, 1e-2) - math.e)
fine = abs(central_difference(D.exponential, 1.0, 5e-3) - math.e)
print("halving_h_quarters_central_error", 3.9 < coarse / fine < 4.1)
coarse_f = abs(forward_difference(D.exponential, 1.0, 1e-2) - math.e)
fine_f = abs(forward_difference(D.exponential, 1.0, 5e-3) - math.e)
print("halving_h_halves_forward_error", 1.9 < coarse_f / fine_f < 2.1)

# -- the U ------------------------------------------------------------------
fe = error_curve(D.exponential, D.U_POINT, D.U_EXACT_SLOPE, D.U_WIDTHS, forward_difference)
ce = error_curve(D.exponential, D.U_POINT, D.U_EXACT_SLOPE, D.U_WIDTHS, central_difference)
print("grid_size", len(D.U_WIDTHS))
print("grid_first", f"{D.U_WIDTHS[0]:.0e}")
print("grid_last", f"{D.U_WIDTHS[-1]:.0e}")
print("forward_curve_is_u", is_u_shaped(fe))
print("central_curve_is_u", is_u_shaped(ce))
bf_h, bf_e = best_step(D.U_WIDTHS, fe)
bc_h, bc_e = best_step(D.U_WIDTHS, ce)
print("best_forward_h", f"{bf_h:.3e}")
print("best_forward_error", f"{bf_e:.6e}")
print("best_central_h", f"{bc_h:.3e}")
print("best_central_error", f"{bc_e:.6e}")
print("best_central_in_band", 1e-7 <= bc_h <= 1e-4)
print("best_forward_in_band", 1e-9 <= bf_h <= 1e-6)
print("best_central_beats_best_forward", bc_e < bf_e)
print("large_h_end_far_worse", ce[0] > 100.0 * bc_e)
print("small_h_end_far_worse", ce[-1] > 100.0 * bc_e)
print("minimum_is_interior", 0 < ce.index(min(ce)) < len(ce) - 1)
print("balance_prediction_forward",
      0.1 < bf_h / math.sqrt(2.0 * D.EPSILON) < 10.0)
print("balance_prediction_central",
      0.1 < bc_h / (3.0 * D.EPSILON) ** (1.0 / 3.0) < 10.0)
print("tiny_h_returns_zero", forward_difference(D.exponential, 1.0, 1e-300))
print("tiny_h_samples_collide", math.exp(1.0 + 1e-300) == math.exp(1.0))
print("epsilon", D.EPSILON == float(np.finfo(np.float64).eps))

# -- stationary points ------------------------------------------------------
H, TOL = D.STATIONARY_WIDTH, D.STATIONARY_TOL
print("parabola_vertex_slope", central_difference(D.parabola, 2.0, H))
print("cubic_min_slope_is_zero", abs(central_difference(D.cubic, 1.0, H)) < TOL)
print("cubic_max_slope_is_zero", abs(central_difference(D.cubic, -1.0, H)) < TOL)
print("cube_step_slope_is_zero", abs(central_difference(D.plain_cube, 0.0, H)) < TOL)
print("parabola_curvature", round(second_difference(D.parabola, 2.0, H), 5))
print("cubic_min_curvature", round(second_difference(D.cubic, 1.0, H), 5))
print("cubic_max_curvature", round(second_difference(D.cubic, -1.0, H), 5))
print("cube_step_curvature", round(second_difference(D.plain_cube, 0.0, H), 5))
print("curvature_separates_min_from_max",
      second_difference(D.cubic, 1.0, H) > 0.0 > second_difference(D.cubic, -1.0, H))
print("classify_parabola", classify_stationary_point(D.parabola, 2.0, H, TOL))
print("classify_cubic_min", classify_stationary_point(D.cubic, 1.0, H, TOL))
print("classify_cubic_max", classify_stationary_point(D.cubic, -1.0, H, TOL))
print("classify_cube_step", classify_stationary_point(D.plain_cube, 0.0, H, TOL))
print("classify_quartic", classify_stationary_point(lambda x: x**4, 0.0, H, TOL))
print("classify_sloping", classify_stationary_point(D.cubic, 0.0, H, TOL))
print("downhill_signs", "|".join(
    "left" if central_difference(D.parabola, x, H) > TOL else "right"
    for x in (-1.0, 0.5, 1.5, 2.5, 4.0)))

# -- corners ----------------------------------------------------------------
print("abs_forward_at_zero", forward_difference(D.absolute, 0.0, D.CORNER_WIDTH))
print("abs_backward_at_zero", backward_difference(D.absolute, 0.0, D.CORNER_WIDTH))
print("abs_central_at_zero", central_difference(D.absolute, 0.0, D.CORNER_WIDTH))
print("abs_central_never_converges",
      all(central_difference(D.absolute, 0.0, h) == 0.0
          for h in (1e-2, 1e-5, 1e-8, 1e-11, 1e-14)))
print("abs_curvature_at_1e3", round(second_difference(D.absolute, 0.0, 1e-3), 1))
print("abs_curvature_at_1e5", round(second_difference(D.absolute, 0.0, 1e-5), 1))
print("relu_forward_at_zero", forward_difference(D.relu, 0.0, D.CORNER_WIDTH))
print("relu_backward_at_zero", backward_difference(D.relu, 0.0, D.CORNER_WIDTH))
print("relu_central_at_zero", central_difference(D.relu, 0.0, D.CORNER_WIDTH))
print("one_sided_gap_at_corner",
      abs(forward_difference(D.absolute, 0.0, D.CORNER_WIDTH)
          - backward_difference(D.absolute, 0.0, D.CORNER_WIDTH)))
print("one_sided_gap_when_smooth",
      abs(forward_difference(D.square, 3.0, D.CORNER_WIDTH)
          - backward_difference(D.square, 3.0, D.CORNER_WIDTH)) < 1e-3)

# -- numpy ------------------------------------------------------------------
print("numpy_gradient_matches_bit_for_bit",
      numpy_gradient_slope(D.exponential, 1.0, D.COMPARE_WIDTH)
      == central_difference(D.exponential, 1.0, D.COMPARE_WIDTH))
print("numpy_coordinates_route_differs",
      numpy_gradient_slope_from_coordinates(D.exponential, 1.0, D.COMPARE_WIDTH)
      != numpy_gradient_slope(D.exponential, 1.0, D.COMPARE_WIDTH))
print("numpy_coordinates_route_differs_slightly",
      abs(numpy_gradient_slope_from_coordinates(D.exponential, 1.0, D.COMPARE_WIDTH)
          - numpy_gradient_slope(D.exponential, 1.0, D.COMPARE_WIDTH)) < 1e-10)
PY
)"

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

check_eq "the car's average speed over six seconds is 24 m/s" "24.0" "$(get car_whole_trip)"
check_eq "over the fourth second alone it is 28 m/s" "28.0" "$(get car_fourth_second)"
check_eq "over the third second it is 20 m/s" "20.0" "$(get car_third_second)"
check_eq "an interval of zero width raises ZeroDivisionError rather than guessing" \
  "ZeroDivisionError" "$(get zero_width_raises)"

check_eq "the shrinking sequence is 7, 6.1, 6.01, 6.001" \
  "7.000000|6.100000|6.010000|6.001000" "$(get settle_sequence)"
check_eq "and every term equals 6 + h to within 1e-12" "True" "$(get settle_matches_six_plus_h)"
check_eq "each interval lands closer to 6 than the one before" "True" "$(get settle_monotone)"
check_eq "approaching from the right comes down from above" "True" "$(get settle_from_above)"
check_eq "approaching from the left comes up from below" "True" "$(get settle_from_below)"
check_eq "the tangent at x = 3 has slope 6" "6.0" "$(get tangent_slope)"
check_eq "and intercept -9, so the line is y = 6x - 9" "-9.0" "$(get tangent_intercept)"

check_eq "all eight derivative rules agree with the arithmetic" "True" "$(get rules_all_within_tolerance)"
check_eq "and their exact values are the eight documented numbers" \
  "0.0000|6.0000|25.3125|-0.2500|30.0000|16.0000|2.7183|0.2500" "$(get rules_exact_values)"
check_eq "the slope of 2**x at zero is the natural log of 2, not 1" \
  "$(get ln_two)" "$(get slope_of_two_to_the_x_at_zero)"
check_eq "the slope of e**x at zero IS 1, which is what makes e special" \
  "1.0" "$(get slope_of_e_to_the_x_at_zero)"

check_eq "the forward difference of x**2 at 3 with h = 0.1 is 6.1" \
  "6.100000000000012" "$(get forward_on_parabola)"
check_eq "the backward difference is 5.9" "5.899999999999999" "$(get backward_on_parabola)"
check_eq "the central difference is exactly 6 on a parabola" "True" "$(get central_on_parabola_is_six)"
check_eq "at h = 1e-5 on e**x the central error is a thousand times smaller" \
  "True" "$(get central_beats_forward_thousandfold)"
check_eq "halving h quarters the central error" "True" "$(get halving_h_quarters_central_error)"
check_eq "halving h only halves the forward error" "True" "$(get halving_h_halves_forward_error)"
echo "  (measured on this run: forward error $(get forward_error_at_1e5), central error $(get central_error_at_1e5), a factor of $(get central_advantage_ratio) -- reported, not asserted)"

check_eq "the error grid holds 27 step sizes" "27" "$(get grid_size)"
check_eq "starting at 1e-1" "1e-01" "$(get grid_first)"
check_eq "and ending at 1e-14" "1e-14" "$(get grid_last)"
check_eq "the forward error curve is U-shaped" "True" "$(get forward_curve_is_u)"
check_eq "the central error curve is U-shaped" "True" "$(get central_curve_is_u)"
check_eq "the minimum is in the interior, not at either end" "True" "$(get minimum_is_interior)"
check_eq "the largest step is more than a hundred times worse than the best" \
  "True" "$(get large_h_end_far_worse)"
check_eq "and so is the smallest, which is the surprising half" \
  "True" "$(get small_h_end_far_worse)"
check_eq "the best central step is in the 1e-7 to 1e-4 band" "True" "$(get best_central_in_band)"
check_eq "the best forward step is in the 1e-9 to 1e-6 band" "True" "$(get best_forward_in_band)"
check_eq "the best central step beats the best forward step" "True" "$(get best_central_beats_best_forward)"
check_eq "the measured forward optimum is within 10x of sqrt(2*EPSILON)" \
  "True" "$(get balance_prediction_forward)"
check_eq "the measured central optimum is within 10x of (3*EPSILON)**(1/3)" \
  "True" "$(get balance_prediction_central)"
check_eq "a step of 1e-300 returns exactly 0.0, silently" "0.0" "$(get tiny_h_returns_zero)"
check_eq "because exp(1 + 1e-300) and exp(1) are the same float64" \
  "True" "$(get tiny_h_samples_collide)"
check_eq "the EPSILON in dataset.py is numpy's float64 epsilon" "True" "$(get epsilon)"
echo "  (measured on this run: best forward h $(get best_forward_h) at error $(get best_forward_error); best central h $(get best_central_h) at error $(get best_central_error) -- reported, not asserted)"

check_eq "the slope at the parabola's vertex is exactly zero" "0.0" "$(get parabola_vertex_slope)"
check_eq "the slope at the cubic's minimum is zero" "True" "$(get cubic_min_slope_is_zero)"
check_eq "the slope at the cubic's maximum is zero too" "True" "$(get cubic_max_slope_is_zero)"
check_eq "and so is the slope at the cubic's flat step" "True" "$(get cube_step_slope_is_zero)"
check_eq "the parabola's curvature is 2" "2.0" "$(get parabola_curvature)"
check_eq "the curvature at the cubic's minimum is +6" "6.0" "$(get cubic_min_curvature)"
check_eq "the curvature at the cubic's maximum is -6" "-6.0" "$(get cubic_max_curvature)"
check_eq "the curvature at the flat step is 0, which decides nothing" \
  "0.0" "$(get cube_step_curvature)"
check_eq "so the second derivative separates the minimum from the maximum" \
  "True" "$(get curvature_separates_min_from_max)"
check_eq "the parabola vertex classifies as a minimum" "minimum" "$(get classify_parabola)"
check_eq "the cubic at +1 classifies as a minimum" "minimum" "$(get classify_cubic_min)"
check_eq "the cubic at -1 classifies as a maximum" "maximum" "$(get classify_cubic_max)"
check_eq "x**3 at 0 classifies as undecided rather than as a minimum" \
  "undecided" "$(get classify_cube_step)"
check_eq "and so does x**4 at 0, which genuinely IS a minimum" \
  "undecided" "$(get classify_quartic)"
check_eq "a point with a real slope is not stationary at all" \
  "not stationary" "$(get classify_sloping)"
check_eq "the sign of the slope points downhill towards the minimum at every x" \
  "right|right|right|left|left" "$(get downhill_signs)"

check_eq "the forward difference of |x| at 0 is +1" "1.0" "$(get abs_forward_at_zero)"
check_eq "the backward difference of |x| at 0 is -1" "-1.0" "$(get abs_backward_at_zero)"
check_eq "the central difference of |x| at 0 is 0.0, where no derivative exists" \
  "0.0" "$(get abs_central_at_zero)"
check_eq "and no smaller h ever reveals a limit that is not there" \
  "True" "$(get abs_central_never_converges)"
check_eq "the curvature at the corner is 2/h: 2,000 at h = 1e-3" "2000.0" "$(get abs_curvature_at_1e3)"
check_eq "and 200,000 at h = 1e-5, so it diverges rather than converging" \
  "200000.0" "$(get abs_curvature_at_1e5)"
check_eq "relu's forward difference at 0 is 1" "1.0" "$(get relu_forward_at_zero)"
check_eq "relu's backward difference at 0 is 0" "0.0" "$(get relu_backward_at_zero)"
# Section 6 re-runs this script with D108_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_relu_central="0.5"
if [ -n "${D108_SELF_TEST:-}" ]; then
  expected_relu_central="1.0"   # the belief that a corner has its right-hand slope
fi
check_eq "relu's central difference at 0 is 0.5, the average of two disagreeing slopes" \
  "${expected_relu_central}" "$(get relu_central_at_zero)"
check_eq "the one-sided rules disagree by 2 at the corner" "2.0" "$(get one_sided_gap_at_corner)"
check_eq "and agree where a derivative really exists" "True" "$(get one_sided_gap_when_smooth)"

check_eq "np.gradient with scalar spacing is our central difference, bit for bit" \
  "True" "$(get numpy_gradient_matches_bit_for_bit)"
check_eq "but passing coordinates instead takes a different arithmetic route" \
  "True" "$(get numpy_coordinates_route_differs)"
check_eq "which differs only in the last few bits" \
  "True" "$(get numpy_coordinates_route_differs_slightly)"

# --------------------------------------------------------------------------
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 -- 1.0, which is what you would believe if you assumed a
# corner simply takes its right-hand slope -- and asserts that the re-run
# reports the failure and exits non-zero. If this section passes, section 5 is
# not decorative.
if [ -z "${D108_SELF_TEST:-}" ]; then
  self_out="$(D108_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: relu's central difference at 0 is 0.5"*)
      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

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 ]