A Python variable is a name tag, not a box. planet = "Mercury" creates a string object somewhere in memory and ties the name planet to it. The object carries the type (int, float, str, bool, None); the name carries nothing — which is why you can rebind the same name to an integer one line later. type(x) asks the object, never the name.

The distinction that quietly explains half of Python's surprising behavior: mutable vs immutable. Numbers, strings, and tuples are immutable — "changing" one actually builds a new object and moves your tag. Lists, dicts, and sets are mutable: they change in place, and every name pointing at that object sees the change. Pass a list into a function and mutate it there — the caller sees it, because the parameter was just another tag on the same box.

• Core scalar types: int, float, str, bool, and None ("no value yet").
• == vs is: value equality vs same-object identity. id(x) reveals the identity.
• Shared mutable arguments are the #1 beginner bug — nothing is copied on assignment.

planet    = "Mercury"      # str  — immutable
gravity_g = 0.38           # float
moons     = 0              # int
habitable = False          # bool
station   = None           # NoneType: "no value yet"

def add_target(missions):
    missions.append("Titan")    # mutates the SHARED object

targets = ["Mars", "Europa"]
backup  = targets               # second tag, same list
add_target(targets)
print(backup)                   # ['Mars', 'Europa', 'Titan'] — both names see it

The list is Python's workhorse container: ordered, mutable, and happy to hold anything. Fluency here is three moves: indexing (neos[0] first, neos[-1] last), slicing (neos[1:4] — start inclusive, stop exclusive; neos[::-1] is a reversed copy), and comprehensions — [n for n in neos if n.hazardous] reads like a sentence and replaces four lines of loop.

Know the performance shape: append is cheap (amortized O(1)), but insert(0, x) shifts every element — O(n). If you need a queue, that's what collections.deque is for.

• Slices clip, never crash: neos[2:999] just returns what exists.
• Slices copy the list, not the elements — top5 = ranked[:5] is a new list of shared references.

neos = ["Apophis", "Bennu", "Ryugu", "Eros", "Itokawa"]

print(neos[0], neos[-1])       # Apophis Itokawa
print(neos[1:4])               # ['Bennu', 'Ryugu', 'Eros']
print(neos[::-1])              # reversed COPY — original untouched

diameters_km = [0.37, 0.49, 0.90, 16.8, 0.33]
big = [n for n, d in zip(neos, diameters_km) if d > 0.4]
print(big)                     # ['Bennu', 'Ryugu', 'Eros']

A tuple is a sealed list: ordered like a list, but immutable after creation. That immutability is not a limitation — it's the feature. Immutable means hashable, and hashable means a tuple can be a dictionary key or live in a set. A list can do neither, because its hash would change every time it mutated.

Tuples are also the hidden mechanism behind returning multiple values: return peri, aph actually returns one tuple, and p, a = orbit() unpacks it. Even the classic swap a, b = b, a is tuple packing and unpacking in disguise.

• The comma makes the tuple, not the parentheses: (3.7,) is a 1-tuple, (3.7) is just a float.
• Immutable ≠ deep-frozen: a tuple holding a list can still have that inner list mutated.

def orbit_extremes(a_au, e):
    peri = a_au * (1 - e)
    aph  = a_au * (1 + e)
    return peri, aph                 # returns ONE tuple

p, a = orbit_extremes(17.8, 0.967)   # Halley's comet — unpacked
print(f"{p:.2f} AU … {a:.2f} AU")     # 0.59 AU … 35.01 AU

stations = {(-35.4, 148.9): "Canberra", (40.4, -4.2): "Madrid"}
print(stations[(-35.4, 148.9)])       # tuples as dict keys — works
# stations[[-35.4, 148.9]] = "?"      # TypeError: unhashable type: 'list'

A dict is a hash map: hash(key) points straight to a slot, so lookup is O(1) average — no scanning. It is the single most important structure for an AI engineer because every JSON payload becomes a dict the moment json.loads runs: NASA API responses, model configs, LLM tool-call arguments — all dicts.

Fluency markers: .get(key, default) instead of KeyError roulette on ragged payloads, .items() to loop key–value pairs, dict comprehensions to reshape data in one line, and the | merge operator (3.9+). Since Python 3.7, dicts preserve insertion order.

• .get("close_approach_data", []) — defaults are how you survive inconsistent API responses.
• Keys must be hashable: str, int, tuple yes; list, dict, set no.

neo = {                              # one object from the NEO feed JSON
    "name": "(2010 PK9)",
    "is_potentially_hazardous_asteroid": True,
    "estimated_diameter": {"kilometers": {"estimated_diameter_max": 0.62}},
}

d_km = neo["estimated_diameter"]["kilometers"]["estimated_diameter_max"]
sentry = neo.get("sentry_data", "not on the Sentry risk list")  # no KeyError

feed  = [{"name": "Bennu", "d_km": 0.49}, {"name": "Ryugu", "d_km": 0.90}]
sizes = {n["name"]: n["d_km"] for n in feed}   # dict comprehension
for name, km in sizes.items():
    print(f"{name}: {km} km")

A set is an unordered bag of unique hashable items with O(1) membership tests — x in my_set hashes; x in my_list scans. Building one is instant deduplication: set(ids) collapses repeats.

The real power is set algebra: & intersection ("in both"), | union ("in either"), - difference ("in A but not B"), ^ symmetric difference ("in exactly one"). Messy dedupe-and-compare logic becomes a one-liner you can read aloud.

• No order: need dedupe that preserves order? list(dict.fromkeys(items)).
• Elements must be hashable — tuples in, lists out (topic 3 pays off again).

survey_a = {"Apophis", "Bennu", "Ryugu", "Eros"}
survey_b = {"Bennu", "Eros", "Itokawa"}

print(survey_a & survey_b)    # both saw:   {'Bennu', 'Eros'}
print(survey_a | survey_b)    # anyone saw: 5 unique objects
print(survey_a - survey_b)    # only A saw: {'Apophis', 'Ryugu'}

ids = ["a1", "b2", "a1", "c3", "b2"]
print(len(set(ids)))          # 3 — dedupe in one call

Python's for iterates any iterable, not just ranges: lists, dicts (over keys), file handles (over lines), generators, API paginators. Idiomatic Python almost never writes for i in range(len(x)) — instead, enumerate() hands you positions and zip() walks two sequences in lockstep.

while handles loops whose end you can't predict: retries, polling a job status, converging computations. Both loop kinds have an else clause that runs only if the loop finished without break — niche, but perfect for "tried everything, all failed" logic.

• enumerate(neos, start=1) — human-friendly numbering without manual counters.
• zip stops at the shortest input — silently. Reach for itertools.zip_longest when that's a bug.

planets = ["Mercury", "Venus", "Earth", "Mars"]
gravity = [0.38, 0.91, 1.00, 0.38]

for i, (name, g) in enumerate(zip(planets, gravity), start=1):
    print(f"{i}. {name}: {g} g")          # lockstep + positions

attempt = 0
while attempt < 3:                        # retry a flaky endpoint
    attempt += 1
    if fetch_donki_flares() is not None:
        break
else:                                     # ran only if no break fired
    print("DONKI still down after 3 tries")

Python functions are first-class objects: assign them to names, store them in dicts, pass them as arguments, return them from other functions. That single fact is the foundation for decorators (topic 12), callback registries, and every agent tool table you'll ever build.

Flexible signatures: *args gathers extra positional arguments into a tuple, **kwargs gathers keyword arguments into a dict — and on the calling side, * and ** unpack. The classic trap: default values are evaluated once, at def time. A mutable default like acc=[] is one shared list across every call — the fix is the None sentinel pattern.

def log_pass(reading, acc=[]):         # BUG: one list, created at def time
    acc.append(reading)
    return acc

print(log_pass("pass 1"))             # ['pass 1']
print(log_pass("pass 2"))             # ['pass 1', 'pass 2']  — surprise!

def log_pass_fixed(reading, acc=None):
    acc = [] if acc is None else acc   # fresh list per call
    acc.append(reading)
    return acc

print(log_pass_fixed("pass 1"))       # ['pass 1']
print(log_pass_fixed("pass 2"))       # ['pass 2']  — independent calls

A class bundles state (attributes) and behavior (methods) into one blueprint. __init__ runs at construction and sets up the instance; self is the instance — Python passes it to every method explicitly, which is why it appears first in every signature.

Dunder methods plug your objects into Python's built-in machinery: __repr__ controls what debugging shows you, __eq__ defines what == means, __len__ powers len(). Once __eq__ exists, two instances holding equal data compare == True while is stays False — value equality vs object identity, the same distinction from topic 1.

class Asteroid:
    def __init__(self, name, d_km):
        self.name = name
        self.d_km = d_km

    def __repr__(self):
        return f"Asteroid({self.name!r}, {self.d_km})"

    def __eq__(self, other):
        return (self.name, self.d_km) == (other.name, other.d_km)

a = Asteroid("Bennu", 0.49)
b = Asteroid("Bennu", 0.49)
print(a)              # Asteroid('Bennu', 0.49) — thanks, __repr__
print(a == b, a is b) # True False — equal value, different objects

Python has no private keyword — encapsulation is by convention: _name means "internal, touch at your own risk" (a leading __name additionally triggers name-mangling). What Python gives you instead of getter/setter boilerplate is @property: callers read and write what looks like a plain attribute, while your validation logic sits invisibly behind it.

The payoff: you can start with a plain attribute today and upgrade it to a guarded property next month without changing a single caller. That's why idiomatic Python skips the Java-style get_x()/set_x() ritual entirely.

class Cryostat:                       # cools an infrared space telescope
    def __init__(self, temp_k):
        self.temp_k = temp_k          # goes through the setter below

    @property
    def temp_k(self):
        return self._temp_k

    @temp_k.setter
    def temp_k(self, value):
        if value < 0:
            raise ValueError("below absolute zero")
        self._temp_k = value

c = Cryostat(6.5)      # JWST's MIRI runs near 6.5 K — fine
c.temp_k = -3          # ValueError — the gate holds

A subclass reuses and extends a parent: it inherits every attribute and method, can override the ones it wants to change, and can reach the parent's version with super() — most commonly super().__init__(…) so the parent still sets up its own state.

When you call a method, Python walks the MRO (method resolution order): instance's class first, then parents in a linearized order you can inspect with ClassName.__mro__. Practical advice: keep hierarchies shallow — one or two levels — and prefer composition when you're tempted to go deeper.

class SpaceProbe:
    def __init__(self, name):
        self.name = name
    def describe(self):
        return f"{self.name}: generic probe"

class Lander(SpaceProbe):
    def __init__(self, name, site):
        super().__init__(name)          # parent still sets up name
        self.site = site
    def describe(self):                  # override
        return f"{self.name}: landed at {self.site}"

print(Lander("Philae", "comet 67P").describe())
print([c.__name__ for c in Lander.__mro__])  # ['Lander', 'SpaceProbe', 'object']

Polymorphism = different classes answering the same method call, each in its own way. Python's flavor is duck typing: the calling code cares that .run() exists, not which class provides it. If it walks like a tool and quacks like a tool, it's a tool — no interface declaration, no isinstance checks.

This is what makes plugin architectures trivially extensible: to add a capability, you write one new class with the right method and drop it in the list. Nothing else changes. (When you want the static checker to enforce the shape, typing.Protocol gives duck typing a spine — that connects to topic 20.)

class NeoTool:
    def run(self, q): return f"asteroids passing near Earth on {q}"

class DonkiTool:
    def run(self, q): return f"solar flares recorded during {q}"

class EonetTool:
    def run(self, q): return f"natural events visible from orbit: {q}"

for tool in [NeoTool(), DonkiTool(), EonetTool()]:
    print(tool.run("2026-08-20"))    # same call, three behaviors

A decorator is a function that takes a function and returns a replacement — usually a wrapper that does something before and/or after calling the original. The @timed line above a def is pure sugar for fetch = timed(fetch): the name now points at the wrapped version.

Decorators are built on closures — the wrapper "remembers" the fn it wraps — and on functions being first-class values (topic 7). Use functools.wraps so the wrapper inherits __name__ and __doc__; without it, logging, debugging, and frameworks all see a function anonymously named wrapper. Every framework annotation you'll ever touch is this mechanism: @app.get, @lru_cache, @retry, @pytest.fixture.

import functools, time

def timed(fn):
    @functools.wraps(fn)               # keep fn's __name__, __doc__
    def wrapper(*args, **kwargs):
        t0 = time.perf_counter()
        result = fn(*args, **kwargs)
        print(f"{fn.__name__}: {time.perf_counter() - t0:.2f}s")
        return result
    return wrapper

@timed
def fetch_neo_feed(days=7):
    time.sleep(0.4)                    # stand-in for the real API call
    return ["Apophis", "Bennu"]

fetch_neo_feed()                       # fetch_neo_feed: 0.40s
print(fetch_neo_feed.__name__)         # 'fetch_neo_feed' — wraps did its job

Put yield in a function and it becomes a generator: calling it runs nothing yet; each request for the next value runs the body to the next yield and freezes there, locals intact. Values are produced one at a time, on demand — constant memory over any size of data.

This is the shape of all streaming: a 10 GB log file line by line, LLM tokens as they're generated, database cursor pages. Generator expressions — (x for x in …) — are the lazy sibling of list comprehensions, and you can chain generators into pipelines where nothing computes until the end consumes.

def thruster_anomalies(path):
    with open(path) as f:
        for line in f:                # file objects are lazy too
            if "THRUSTER" in line:
                yield line.rstrip()   # freeze here until next request

events = thruster_anomalies("voyager_telemetry.log")
first = next(events)                  # nothing was loaded up front
for line in events:                   # resumes exactly where it froze
    print(line)

lengths = (len(l) for l in events)   # generator expression: still lazy

The iterator protocol is the machinery under every for loop: iter(x) calls x.__iter__() to get an iterator, then the loop calls __next__() repeatedly until it raises StopIteration. A for loop is sugar for exactly that while/try dance — write one by hand once and the magic becomes mundane.

Key vocabulary: an iterable has __iter__ (a list — you can loop it many times); an iterator also has __next__ and is consumed once. Generators (topic 13) are iterators the language writes for you — which is why looping a generator a second time yields nothing.

class Countdown:
    def __init__(self, start):
        self.n = start
    def __iter__(self):
        return self                    # I am my own iterator
    def __next__(self):
        if self.n == 0:
            raise StopIteration        # how every loop actually ends
        self.n -= 1
        return self.n + 1

for t in Countdown(3):
    print(f"T-minus {t}")              # T-minus 3, 2, 1 … liftoff

The with statement guarantees setup and teardown even when exceptions fly: files close, locks release, connections return to the pool. Under the hood, __enter__ runs before the block and __exit__ always runs after — it even receives the exception info and can choose to suppress it.

The quick way to write your own: @contextlib.contextmanager on a generator — everything before yield is the enter, everything after (inside a finally) is the exit. One yield, and you have a bulletproof resource guard.

from contextlib import contextmanager
import time

@contextmanager
def stage(name):
    t0 = time.perf_counter()
    try:
        yield                          # the with-block body runs here
    finally:
        print(f"{name}: {time.perf_counter() - t0:.2f}s")

with stage("embed batch"):
    embed(chunks)                      # timed even if embed() raises

The full shape is try/except/else/finally: except catches (name specific exception classes — a bare except: swallows bugs you wanted to see), else runs only when nothing was raised, finally runs no matter what. Catch what you can handle; let everything else travel up.

Two production habits: custom exception classes make failure modes explicit — DonkiTimeout and DonkiBadPayload deserve different handling — and raise … from err preserves the causal chain, so the traceback shows the original failure underneath your domain error. Production LLM code is mostly error handling: timeouts, rate limits, and malformed JSON from the model itself.

import time

class DonkiError(Exception):
    """Space-weather data could not be fetched."""

def fetch_flares_with_retry(day, tries=3):
    for attempt in range(1, tries + 1):
        try:
            return fetch_donki(day)              # the flaky call
        except TimeoutError as err:
            if attempt == tries:
                raise DonkiError(f"gave up after {tries} tries") from err
            time.sleep(2 ** attempt)             # backoff: 2s, 4s

Threads live in one process and share memory — cheap to start, trivial to share data between. But Python's GIL (global interpreter lock) means only one thread executes Python bytecode at any instant, so threads cannot parallelize computation.

The nuance that makes threads useful anyway: a thread that is waiting on IO — a socket, a disk read, a sleep — releases the GIL. While one thread waits on NASA's servers, another runs. So for IO-bound work, ten threads ≈ ten requests in flight and a near-10× wall-clock win; for CPU-bound work, threads add overhead and win nothing.

The practical interface is concurrent.futures.ThreadPoolExecutor: pool.map for uniform fan-out, submit + as_completed when you want results as they land. Shared memory cuts both ways — two threads mutating one dict is a race. The sane pattern: threads share nothing, each returns a value, and the pool collects them.

IO-bound vs CPU-bound is the one-question litmus for all three concurrency topics: is the CPU waiting, or working? Waiting on network/disk → threads or asyncio. Working through math/parsing → processes (topic 18). Answer this before choosing a tool, every time.

from concurrent.futures import ThreadPoolExecutor
import requests, time

DAYS = [f"2026-08-{d:02d}" for d in range(1, 11)]

def fetch(day):                        # pure IO: the GIL releases while waiting
    r = requests.get(NEO_URL, params={"start_date": day, "api_key": KEY})
    return day, r.json()["element_count"]

t0 = time.perf_counter()
with ThreadPoolExecutor(max_workers=10) as pool:
    counts = list(pool.map(fetch, DAYS))   # 10 requests in flight at once
print(f"{time.perf_counter() - t0:.1f}s")  # ≈ one round-trip, not ten

Separate processes mean separate interpreters, separate memory, and — crucially — separate GILs. This is how Python gets true parallelism across CPU cores, and it's the tool for CPU-bound work: feature extraction, parsing, number crunching.

The costs are real: each process must start up (fork/spawn), and every argument and result crosses the boundary by pickling — serialized, shipped, deserialized. Ship a huge NumPy array to a worker for a tiny computation and the shipping eats the win. The rule: the compute must dwarf the transfer. That's why Pool.map over ten million numbers wins big, while the same call on a ten-element list loses to a plain loop.

• multiprocessing.Pool.map (or ProcessPoolExecutor — same API as the thread pool, swap one word).
• Send paths, not arrays: let workers load their own data to dodge the pickling tax.
• Guard the entry point with if __name__ == "__main__": — on spawn platforms (macOS, Windows) workers re-import your module.

from multiprocessing import Pool
import numpy as np

def light_curve_features(path):        # CPU-bound: FFT + stats per star
    flux = np.load(path)
    return path, np.abs(np.fft.rfft(flux)).max(), flux.std()

if __name__ == "__main__":             # required on macOS/Windows spawn
    paths = [f"curves/kic_{i}.npy" for i in range(100_000)]
    with Pool(processes=8) as pool:
        feats = pool.map(light_curve_features, paths, chunksize=64)
    # paths cross the boundary (tiny); arrays are loaded inside workers

Asyncio is single-threaded cooperative concurrency. An async def function is a coroutine: it runs until it hits await on something that isn't ready (a network response, a timer), then hands control back to the event loop, which immediately runs whatever else is ready. When the IO completes, the coroutine resumes right after its await. No threads, no GIL wrestling — one thread juggling thousands of open requests.

This is THE core skill for LLM engineering, because LLM systems are IO-bound end to end: model API calls, vector DB queries, tool HTTP calls. The two verbs to master: asyncio.gather(*coros) fans out many awaits and collects all results; asyncio.Semaphore(n) caps how many run at once — backpressure, so 500 concurrent calls don't melt a rate limit. asyncio.wait_for adds timeouts.

The discipline: inside a coroutine, never call blocking code. time.sleep, requests.get, or a heavy computation freezes the entire loop — every other task stalls. Use asyncio.sleep, httpx.AsyncClient, and push CPU work to an executor. FastAPI endpoints are coroutines running on exactly this loop: one blocking call in one endpoint degrades every request in flight.

• await means: "park me here; loop, go run someone else."
• gather = fan-out · Semaphore = backpressure · wait_for = timeout.
• asyncio.run(main()) starts the loop; you only write coroutines.

Threads vs asyncio, both IO-bound tools: threads are preemptive (the OS switches them anywhere, races possible, ~MB of stack each); coroutines are cooperative (switches happen only at visible awaits, ~KB each — thousands are cheap). Choose asyncio when you control the code and it's IO; use threads when you're stuck with a blocking library — or bridge with asyncio.to_thread.

import asyncio, httpx

SEM = asyncio.Semaphore(5)                  # backpressure: max 5 in flight

async def fetch(client, day):
    async with SEM:                         # waits here if 5 are already out
        r = await client.get(NEO_URL, params={"start_date": day})
        return day, r.json()["element_count"]

async def main():
    days = [f"2026-08-{d:02d}" for d in range(1, 21)]
    async with httpx.AsyncClient(timeout=10) as client:
        counts = await asyncio.gather(*(fetch(client, d) for d in days))
    print(counts)                           # 20 calls ≈ 4 round-trips, 1 thread

asyncio.run(main())

Type hints are annotations — str, list[int], dict[str, float], str | None (the modern Optional), TypedDict for dict shapes, Callable, generics — that document intent and let mypy/pyright catch bugs before anything runs. The interpreter mostly ignores them; the tooling and your teammates do not.

In modern AI codebases hints are load-bearing, not optional: FastAPI turns your parameter annotations into request parsing, validation, and OpenAPI docs; Pydantic turns them into runtime validation of LLM and API payloads. Reading a well-typed signature tells you the whole contract without opening the body.

• Start at the boundaries: public functions, API payloads, tool schemas — where wrong types cost the most.
• X | None forces the question every caller must answer: what happens when it's None?

from typing import TypedDict

class Neo(TypedDict):
    name: str
    d_km: float
    hazardous: bool

def riskiest(neos: list[Neo], min_d: float = 0.3) -> Neo | None:
    big = [n for n in neos if n["d_km"] >= min_d and n["hazardous"]]
    return max(big, key=lambda n: n["d_km"]) if big else None

# $ mypy asteroids.py
# error: Argument "min_d" to "riskiest" has incompatible type "str";
#        expected "float"        ← caught before runtime, for free

A virtual environment is a per-project interpreter plus its own site-packages. Project A keeps numpy 1.26, project B keeps numpy 2.1, and neither knows the other exists. "Activating" just points your shell's PATH at the project's python and pip.

python -m venv .venv is the standard; uv is the fast modern replacement worth adopting now — one tool covering venv + pip + lockfiles, with a resolver measured in milliseconds. The iron rule either way: never install into the system Python.

$ python -m venv .venv
$ source .venv/bin/activate
(.venv) $ pip install httpx==0.27.0
(.venv) $ pip freeze | head -1        # httpx==0.27.0 — lives here only
(.venv) $ deactivate
$ python -c "import httpx"           # ModuleNotFoundError — isolation works

# the modern way — same idea, dramatically faster:
$ uv venv && uv pip install httpx

pyproject.toml is how Python code becomes an installable unit: name, version, dependencies, and build backend in one declarative file. Once a project has one, pip install -e . gives you an editable install — the package is importable from anywhere on the machine, and your edits take effect immediately, no reinstall.

This kills the sys.path.append("../../src") hack for good. Bonus structure: optional dependency groups ([project.optional-dependencies]) keep dev tools like pytest and mypy out of production installs, and a src/ layout keeps accidental imports honest.

# pyproject.toml
[project]
name = "orbital-tools"
version = "0.3.0"
requires-python = ">=3.11"
dependencies = ["httpx>=0.27", "pydantic>=2.7"]

[project.optional-dependencies]
dev = ["pytest", "mypy", "ruff"]

[build-system]
requires = ["hatchling"]
build-backend = "hatchling.build"

# $ pip install -e ".[dev]"
# → 'import orbital_tools' now works from any directory, edits live

The universal skeleton is arrange–act–assert: set up inputs, run the code under test, assert on the result. Good tests are isolated (no shared state between tests), deterministic (no live network, no randomness without seeds), and named after the behavior they pin down — a failing name should read like a bug report.

Test the contract, not the implementation — and always cover the edges: zero, empty, and invalid input that should raise (that's what assertRaises is for). This habit is what separates engineers from notebook users.

import unittest
from hazard import hazard_score

class TestHazardScore(unittest.TestCase):
    def test_big_close_asteroid_scores_high(self):
        self.assertGreater(hazard_score(d_km=1.2, miss_ld=1.5), 0.8)

    def test_zero_diameter_scores_zero(self):
        self.assertEqual(hazard_score(d_km=0.0, miss_ld=5.0), 0.0)

    def test_negative_distance_raises(self):
        with self.assertRaises(ValueError):
            hazard_score(d_km=0.4, miss_ld=-1.0)

# $ python -m unittest — normal case, edge case, error case

Pytest is the de-facto runner: plain assert statements (it introspects failures and shows you both sides), auto-discovery of test_*.py, and two superpowers. Fixtures declare shared setup once — any test that names the fixture as a parameter receives it, composed and scoped (function/module/session). @pytest.mark.parametrize runs one test body over many input/expected pairs, each reported as its own case.

Built-in fixtures like monkeypatch (patch env vars and attributes), capsys, and tmp_path cover most mocking needs. Every serious Python AI repo you'll touch runs pytest.

import pytest
from hazard import hazard_score

@pytest.fixture
def apophis():
    return {"d_km": 0.37, "miss_ld": 0.098}   # 2029 close approach

def test_apophis_is_flagged(apophis):          # fixture injected by name
    assert hazard_score(**apophis) > 0.5

@pytest.mark.parametrize("d_km,miss_ld,expected", [
    (0.0,  5.0,  0.0),
    (1.2,  1.5,  0.97),
    (0.05, 60.0, 0.01),
])
def test_known_scores(d_km, miss_ld, expected):
    assert hazard_score(d_km, miss_ld) == pytest.approx(expected, abs=0.02)

A NumPy ndarray is a typed, contiguous block of memory, and operations on it run as compiled C loops — routinely 10–100× faster than Python loops. The mental shift is to think in whole arrays: d_km / miss divides every pair at once; no index in sight.

The concept to actually master is broadcasting — how mismatched shapes stretch to fit. Align shapes from the right: each dimension must be equal, or 1 (it stretches), or missing (it's added). So (3,1) + (1,4) yields (3,4). Master this here, because PyTorch tensors follow exactly the same rules — every shape error you'll ever debug in deep learning is this topic.

import numpy as np

d_km = np.array([0.37, 0.49, 0.90, 16.8])    # NEO diameters
miss = np.array([0.098, 0.48, 6.1, 42.0])    # miss distance (lunar dist.)

danger = d_km / miss                          # vectorized: every pair at once
print(danger.round(2))                        # [3.78 1.02 0.15 0.4 ]

col = np.arange(3).reshape(3, 1)             # shape (3,1)
row = np.arange(4).reshape(1, 4) * 10        # shape (1,4)
print((col + row).shape)                      # (3,4) — predicted before running?

A DataFrame is a table with labeled columns, and pandas is how tabular data gets interrogated: boolean-mask filtering (df[df["d_km"] > 0.5]), groupby — the split-apply-combine pattern: split rows into groups, apply an aggregate, combine into a small result — merge (SQL joins between frames), and missing-value handling (.isna, .fillna, .dropna).

Every dataset you'll ever evaluate a model on passes through pandas first — column operations are vectorized NumPy underneath, so the same "whole columns, not loops" instinct applies.

import pandas as pd

df = pd.read_csv("neo_feed.csv")              # one row per close approach
big = df[df["d_km"] > 0.5]                    # boolean-mask filter

per_day = (df.groupby("date")["d_km"]
             .agg(count="size", max_d="max")   # split → apply
             .reset_index())                    # → combine

weather = pd.read_csv("donki_flares.csv")      # second table
full = per_day.merge(weather, on="date", how="left")
full["flare_class"] = full["flare_class"].fillna("quiet")

Matplotlib's real API is the object model: a Figure is the canvas; each Axes is one plot living on it. fig, ax = plt.subplots() hands you both, and from there everything is explicit: ax.plot, ax.set_xlabel, ax.legend. The plt.plot(…) state machine is fine for a throwaway — the Axes API is how you control subplots, labels, and legends when it matters.

plt.subplots(2, 2) gives you a grid of Axes for side-by-side comparisons, and fig.savefig("out.png", dpi=150) makes the artifact. Every training curve you will ever debug lives on an Axes.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(7, 4))          # the object model
ax.plot(epochs, train_loss, label="train")
ax.plot(epochs, val_loss, label="val", linestyle="--")
ax.scatter([best_epoch], [min(val_loss)], color="red", zorder=3)
ax.set_xlabel("epoch")
ax.set_ylabel("loss")
ax.set_title("Exoplanet CNN — learning curve")
ax.legend()
fig.savefig("learning_curve.png", dpi=150)

Plotly renders charts as interactive HTML/JS: hover for per-point detail, zoom and pan into dense regions, click legend entries to toggle series. plotly.express gives you one-liner figures from a DataFrame; graph_objects is the lower-level control layer when you outgrow it.

The tool choice is about the audience: a human exploring results wants Plotly's interactivity; a paper, README, or CI artifact wants matplotlib's static image. fig.write_html("chart.html") produces a self-contained file anyone can open — no Python required on their end.

import plotly.express as px

fig = px.scatter(
    df, x="miss_distance_ld", y="d_km",
    color="hazardous", log_x=True,
    hover_name="name",                          # big label on hover
    hover_data={"velocity_kps": ":.1f", "date": True},
    title="This week's close approaches",
)
fig.write_html("neo_week.html")                # self-contained, shareable

Seaborn is statistical plotting sugar over matplotlib: heatmap (correlation matrices with annotations), pairplot (every feature against every other), histplot/kdeplot (distributions), boxplot (group comparisons) — each one line from a DataFrame, with sensible statistical defaults.

Because every seaborn function draws onto a matplotlib Axes and returns it, everything from topic 27 still applies: pass ax=, tweak titles and limits, save with fig.savefig. It's the speed tool for EDA — exploration before modeling.

import seaborn as sns
import matplotlib.pyplot as plt

num = df[["transit_depth", "duration_h", "stellar_radius", "period_days"]]

fig, ax = plt.subplots(figsize=(6, 5))
sns.heatmap(num.corr(), annot=True, fmt=".2f",
            cmap="vlag", center=0, ax=ax)      # correlations, annotated
ax.set_title("Exoplanet feature correlations")
fig.savefig("corr.png", dpi=150)

Scikit-learn is the classical-ML library, but the durable lesson is its API contract: every estimator learns with .fit(X, y), predicts with .predict(X), and every preprocessor transforms with .transform(X). Learn the rhythm once and you can operate any model in the library — and most tools outside it, because XGBoost, LightGBM, and countless wrappers imitate this exact interface.

Pipeline chains preprocessing and model into one estimator, which quietly fixes the classic leakage bug: the scaler is fitted only on training data, then applied to test data — automatically, every time, including inside cross-validation.

from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression

X_tr, X_te, y_tr, y_te = train_test_split(     # X: transit features
    X, y, test_size=0.2, stratify=y, random_state=42)

clf = Pipeline([
    ("scale", StandardScaler()),               # fitted on X_tr ONLY
    ("model", LogisticRegression(max_iter=500)),
])
clf.fit(X_tr, y_tr)                             # the whole chain learns
print(f"accuracy: {clf.score(X_te, y_te):.3f}") # held-out evaluation

A vector is an ordered list of numbers representing a point or a direction in space. In AI, everything becomes a vector: a planet's features [mass, radius, orbital_period], a sentence's meaning, a user's behavior. Embeddings — the currency of modern AI — are just vectors, and "similar meaning" becomes "small angle between vectors." Once data is a vector, geometry becomes your toolkit: distance, angle, projection.

• Magnitude: the vector's length — np.linalg.norm(v). Often normalized away so only direction matters.
• Direction: where it points — for embeddings, this is what carries the meaning.
• Dimensions: a 3-feature planet lives in 3-D; an OpenAI embedding lives in 1536-D. Same math, more axes.

import numpy as np

flare = np.array([0.8, 0.3, 0.1])   # embedding of "solar flare"
cme   = np.array([0.7, 0.4, 0.2])   # embedding of "CME event"

cos = flare @ cme / (np.linalg.norm(flare) * np.linalg.norm(cme))
print(f"similarity: {cos:.3f}")   # close to 1 → related concepts

A matrix is a grid of numbers that acts as a transformation: multiply a vector by a matrix and the vector gets rotated, stretched, sheared, or projected into a different space. That's the real mental model — not "a table of numbers" but "a machine that moves vectors." Every neural-network layer's weights are a matrix doing exactly this: the layer takes your 4-feature planet vector and transforms it into a new representation where the next layer's job is easier.

• An m×n matrix maps n-dimensional vectors to m-dimensional ones — it can change dimensionality.
• The columns of a matrix tell you where the basis vectors land — the whole transformation follows from that.
• Rotation, scaling, projection: all just special patterns of numbers in the grid.

import numpy as np

theta = np.deg2rad(45)
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])   # 2×2 rotation matrix

probe = np.array([1.0, 0.0])    # probe position on the x-axis
print(R @ probe)                # [0.707 0.707] → rotated 45°

Rows of the left times columns of the right: each cell of the result is one dot product between a row of A and a column of B. An (m×n) matrix times an (n×p) matrix gives an (m×p) result — and the inner dimensions must match, because a row of length n can only be dotted with a column of length n. This one operation is what GPUs spend their lives doing: a transformer forward pass is a long chain of matrix multiplications, and "FLOPs" mostly counts these.

• Shape rule: (m,n) @ (n,p) → (m,p) — read it left to right, the n's cancel.
• Not commutative: A @ B ≠ B @ A in general — order is meaning.
• Cost is O(m·n·p) multiplications — why model size × sequence length dominates your GPU bill.

import numpy as np

planets = np.random.rand(3, 4)   # 3 planets × 4 features   (m×n)
weights = np.random.rand(4, 2)   # 4 features → 2 hidden units (n×p)

hidden = planets @ weights       # (3×4)·(4×2) → (3×2)
print(hidden.shape)              # (3, 2)

bad = np.random.rand(5, 2)
planets @ bad                    # ValueError: shapes (3,4) and (5,2)
                                 # — the inner 4 and 5 don't match

Multiply element-wise and sum: a·b = a₁b₁ + a₂b₂ + … — one number measuring how aligned two vectors are. Positive means pointing the same general way, zero means orthogonal (no relationship along these axes), negative means opposed. This is the raw score inside both attention (query · key = "how much should this token attend to that one") and vector search (query · document = relevance).

• Geometric identity: a·b = |a||b|cos(θ) — alignment scaled by lengths.
• In NumPy: a @ b or np.dot(a, b).

Dot product vs cosine similarity: cosine is just the dot product after normalizing both vectors to length 1 — pure angle, lengths ignored. Vector DBs let you pick either; if your embeddings are already normalized (many are), they rank identically.

import numpy as np

v_sat   = np.array([7.5, 0.0])    # satellite velocity, km/s
v_along = np.array([3.0, 0.0])    # thrust along the orbit
v_ortho = np.array([0.0, 3.0])    # thrust perpendicular
v_retro = np.array([-3.0, 0.0])   # retrograde burn

for name, u in [("along", v_along), ("ortho", v_ortho), ("retro", v_retro)]:
    print(name, v_sat @ u)         # 22.5 · 0.0 · -22.5

Most vectors get knocked off their direction when a matrix transforms them. But for special directions, the matrix acts as pure scaling: Av = λv — the output is the same vector, just stretched or shrunk. The eigenvalue λ is that scaling factor, and the set of eigenvalues summarizes what a transformation fundamentally does: which directions it amplifies (|λ| > 1), which it dampens (|λ| < 1), which it flips (λ < 0).

• Read Av = λv as: "along direction v, the matrix A is just multiplication by the number λ."
• Big eigenvalue spread (condition number) = a hard optimization landscape — connects to topic 43.

import numpy as np

A = np.array([[3.0, 1.0],
              [0.0, 2.0]])
vals, vecs = np.linalg.eig(A)

for lam, v in zip(vals, vecs.T):
    print(np.allclose(A @ v, lam * v))   # True, True — Av = λv holds
print(vals)                              # [3. 2.] — the pure scale factors

The eigenvectors are the special directions from the previous rung — the vectors whose direction survives the transformation (only their length changes, by factor λ). Where eigenvalues tell you "how much", eigenvectors tell you "along which axes". The flagship application: PCA's principal components are exactly the eigenvectors of the data's covariance matrix — the directions along which your data varies most, ranked by their eigenvalues.

• Eigenvector 1 (largest λ): the single direction that captures the most variance.
• Eigenvectors of a symmetric matrix (like covariance) are orthogonal — they form a new axis system for your data.

import numpy as np

# correlated orbital measurements: 500 samples, 2 features
X = np.random.randn(500, 2) @ np.array([[3.0, 0.0], [1.0, 0.5]])
X -= X.mean(axis=0)                 # center the cloud

cov = np.cov(X.T)                    # 2×2 covariance matrix
vals, vecs = np.linalg.eig(cov)
pc1 = vecs[:, np.argmax(vals)]       # direction of most variance
print(pc1)                           # PCA's first principal component

Any matrix — square or not — factors into rotation × scaling × rotation: A = UΣVᵀ. U and V are rotations (orthogonal matrices); Σ is diagonal, holding the singular values in decreasing order — how strongly A stretches along each of its natural axes. The killer property: keeping only the top-k singular values gives the best possible rank-k approximation of A. That is the exact intuition behind LoRA: a weight update can be well-approximated by a low-rank matrix, so you train two skinny matrices instead of the full grid.

• Eigendecomposition needs a square matrix; SVD works on anything — it's the workhorse.
• Singular values decay fast in real data → most matrices are "secretly low-rank".

import numpy as np

T = np.random.randn(200, 64)      # 200 timesteps × 64 telemetry channels
U, S, Vt = np.linalg.svd(T, full_matrices=False)

k = 5                              # keep only the 5 strongest directions
T_k = U[:, :k] * S[:k] @ Vt[:k, :] # best rank-5 approximation of T

err = np.linalg.norm(T - T_k) / np.linalg.norm(T)
print(f"rank-{k} reconstruction error: {err:.1%}")

The derivative is the instantaneous rate of change — the slope of a function at a point. Zoom in far enough on any smooth curve and it looks like a straight line; the derivative is that line's slope. This is the single most important idea in this volume, because all of training is "compute the slope of the loss, walk downhill". If the derivative of loss with respect to a weight is positive, decreasing the weight decreases the loss — so you decrease it.

• Notation: f′(x) or df/dx — "how fast does f change if I nudge x?"
• Numerical check: (f(x+h) − f(x)) / h for tiny h — how autograd is unit-tested.
• f′(x) = 0 marks flat spots: minima, maxima, saddle points — where training wants to land.

def f(x):
    return x ** 2          # e.g. atmospheric drag vs probe velocity

x = 3.0
for h in [1e-1, 1e-3, 1e-6]:
    slope = (f(x + h) - f(x)) / h
    print(f"h={h:.0e}  slope ≈ {slope:.6f}")   # converges on 6.0 = 2x

When a function has several inputs, a partial derivative is the slope with respect to one variable while all others hold still: ∂f/∂x asks "if I nudge only x, how does f respond?" Loss functions in ML take millions or billions of parameters as input; training asks for the loss's partial derivative with respect to every single one — each weight gets its own private answer to "which way should I move?"

• Symbol ∂ ("del") instead of d signals: other variables are frozen.
• Computing one partial = ordinary 1-D differentiation with the other variables treated as constants.

def loss(w, b):                      # tiny model: radius = w·mass + b
    pred = w * 2.0 + b               # one planet, mass = 2.0
    return (pred - 1.4) ** 2         # true radius 1.4

h = 1e-6
dw = (loss(3 + h, 1) - loss(3, 1)) / h   # ∂loss/∂w — b held still
db = (loss(3, 1 + h) - loss(3, 1)) / h   # ∂loss/∂b — w held still
print(dw, db)                              # 22.4, 11.2 — each knob's own slope

The derivative of nested functions multiplies through the layers: d/dx f(g(x)) = f′(g(x)) · g′(x). Rates of change compose by multiplication — if u changes 2× as fast as x, and y changes 3× as fast as u, then y changes 6× as fast as x. A neural network is nothing but deeply nested functions: loss(softmax(W₃·relu(W₂·relu(W₁·x)))). This one rule, applied backward through the network layer by layer, IS backpropagation.

• Each layer contributes one local derivative factor; the gradient at any depth is the product of everything above it.
• Long chains of small factors → vanishing gradients; of large factors → exploding gradients. Both are chain-rule arithmetic.

import numpy as np

def f(x):
    return np.sin(x ** 2)            # nested: sin ∘ square

def f_prime(x):                       # chain rule, by hand
    return np.cos(x ** 2) * 2 * x     # outer′(inner) · inner′

x, h = 1.3, 1e-7
numeric = (f(x + h) - f(x)) / h
print(np.isclose(f_prime(x), numeric))   # True — the rule holds

The gradient ∇f stacks all the partial derivatives into one vector — and that vector points in the direction of steepest increase of the function. So training steps in the opposite direction: w ← w − η·∇loss, where η is the learning rate. "Gradient descent" is literally just that sentence, repeated until the loss stops improving.

• Gradient = vector of partials: ∇f = [∂f/∂x, ∂f/∂y, …] — direction AND steepness (its norm).
• η too small: crawling. η too large: overshooting the valley and diverging. The classic first hyperparameter.
• SGD = same walk, but each step estimates ∇ from a mini-batch instead of the full dataset.

import numpy as np

pos = np.array([3.0, 4.0])       # start high on the loss bowl
lr = 0.1
for step in range(25):
    grad = 2 * pos               # ∇f for f(x,y) = x² + y²
    pos -= lr * grad             # walk downhill

print(pos)                       # ≈ [0.011 0.015] — near the minimum

When a function maps vectors to vectors (n inputs → m outputs), one number can't describe its slope. The Jacobian is the m×n matrix of every output's partial derivative with respect to every input — row i answers "how does output i respond to each input?" It's the bookkeeping structure that makes backprop through whole layers coherent: each layer's local derivative is a Jacobian, and the chain rule becomes Jacobian matrix multiplication.

• Shape: outputs × inputs. A layer mapping 64 → 32 features has a 32×64 Jacobian.
• The Jacobian is also the best linear approximation of the function near a point — "locally, every smooth map is a matrix."

import numpy as np

def orbit_xy(r, theta):               # polar orbit coords → cartesian
    return np.array([r * np.cos(theta), r * np.sin(theta)])

r, th, h = 7000.0, 0.6, 1e-6          # LEO radius (km), true anomaly
J = np.column_stack([
    (orbit_xy(r + h, th) - orbit_xy(r, th)) / h,   # column 1: ∂(x,y)/∂r
    (orbit_xy(r, th + h) - orbit_xy(r, th)) / h,   # column 2: ∂(x,y)/∂θ
])
print(J)   # [[cosθ, -r·sinθ], [sinθ, r·cosθ]] — matches the analytic Jacobian

The Hessian is the matrix of second derivatives — the curvature of the loss surface. Where the gradient says "which way is down", the Hessian says what shape the terrain is: a bowl (all eigenvalues positive — a minimum), a dome (all negative — a maximum), or a saddle point (mixed signs — down in one direction, up in another). It also explains why optimizers like Adam adapt step sizes per direction: when curvature is steep one way and shallow another, a single global learning rate is always wrong for one of them.

• Shape: n×n for n parameters — why you never materialize it for a billion-parameter model.
• Condition number = largest/smallest eigenvalue; big ratio = elongated valley = slow, zig-zagging SGD.

import numpy as np

# f(x, y) = x² + 3y² — an elongated loss valley
H = np.array([[2.0, 0.0],    # ∂²f/∂x² = 2
              [0.0, 6.0]])   # ∂²f/∂y² = 6 → 3× steeper in y

vals = np.linalg.eigvalsh(H)
print(vals)              # [2. 6.] — curvature per principal direction
print(vals[1] / vals[0]) # 3.0 condition ratio: one global lr can't fit both

P(A|B) = P(B|A) · P(A) / P(B): flip a conditional probability around using prior knowledge. It's the mental model for updating beliefs with evidence: start with a prior (how common is A at all?), observe evidence B, and get the posterior (how likely is A now?). Its most famous lesson is the base-rate fallacy: when A is rare, even a very accurate test produces mostly false positives — the tiny prior dominates.

• Prior P(A): belief before evidence. Likelihood P(B|A): how well the evidence fits. Posterior P(A|B): belief after.
• Naive Bayes classifiers, spam filters, and medical-test intuition are all this one equation.

# transit survey: 1% of stars host a detectable planet
p_planet = 0.01
p_signal_given_planet = 0.95    # detector sensitivity
p_signal_given_none   = 0.05    # false-alarm rate

p_signal = (p_signal_given_planet * p_planet
            + p_signal_given_none * (1 - p_planet))
posterior = p_signal_given_planet * p_planet / p_signal
print(f"P(planet | signal) = {posterior:.1%}")   # 16.1% — not 95%!

P(A|B) is the probability of A once you know B happened — you've shrunk the universe to only the B cases and re-measured A inside it: P(A|B) = P(A∩B) / P(B). This "shrink the universe, renormalize" move is the heart of probabilistic reasoning, and everything in probabilistic ML is conditional on something: a language model computes P(next token | all previous tokens); your classifier computes P(planet | features).

• Conditioning ≠ causation: P(A|B) can be high because B causes A, A causes B, or something causes both.
• If P(A|B) = P(A), knowing B told you nothing — that's the definition of independence.

import numpy as np
rng = np.random.default_rng(0)

d1 = rng.integers(1, 7, 100_000)     # two dice, 100k rolls
d2 = rng.integers(1, 7, 100_000)

given = d1 == 3                       # shrink the universe: first die = 3
p = np.mean(d1[given] + d2[given] == 8)   # P(sum = 8 | d1 = 3)
print(f"{p:.4f}  vs analytic 1/6 = {1/6:.4f}")

A random variable is a variable whose value is the numeric outcome of a random process — it's the bridge between real-world randomness and math you can compute with. Flip 10 coins: "the number of heads" is a random variable X taking values 0–10. Random variables come in two flavors: discrete (countable outcomes — dice, event counts) and continuous (any value in a range — heights, temperatures, noise voltages).

• The random variable is the mapping from outcomes to numbers — not the number itself.
• Once something is a random variable, it earns a distribution, a mean, a variance — the next three topics.

import numpy as np
rng = np.random.default_rng(1)

# X = flares detected in a 10-day window, p = 0.5 per day
detections = rng.random((10_000, 10)) < 0.5
X = detections.sum(axis=1)          # 10,000 samples of the random variable

vals, counts = np.unique(X, return_counts=True)
print(dict(zip(vals, counts / 10_000)))   # X's empirical distribution

A distribution is the full description of a random variable: which values are how likely. For discrete variables it's a pmf (probability mass function — a bar per value, bars sum to 1); for continuous ones a pdf (probability density function — a curve where area, not height, is probability); and the cdf (cumulative) answers "what fraction falls at or below x?" This vocabulary is everywhere in AI: a softmax output is a pmf over classes, sampling temperature reshapes that pmf, and a language model is a distribution over next tokens.

import numpy as np
rng = np.random.default_rng(2)

uniform  = rng.uniform(0, 1, 10_000)        # sensor jitter: flat pdf
normal   = rng.normal(0, 1, 10_000)         # thermal noise: bell pdf
binomial = rng.binomial(10, 0.3, 10_000)    # detections per 10 passes: pmf

for name, x in [("uniform", uniform), ("normal", normal), ("binomial", binomial)]:
    print(f"{name:9s} mean={x.mean():.2f}  var={x.var():.2f}")

The bell curve, defined entirely by two numbers: mean μ (where it's centered) and standard deviation σ (how spread out). Its superpower is the Central Limit Theorem: sum (or average) many small independent effects — whatever their individual distributions — and the result tends toward a Gaussian. That's why it shows up everywhere: measurement noise, sensor error, aggregated returns. And it's why weight initialization and noise models default to it in deep learning.

• Notation: X ~ N(μ, σ²). Standard normal: μ=0, σ=1.
• CLT is the reason "average of many measurements" is such a trustworthy quantity.

import numpy as np
rng = np.random.default_rng(3)

# CLT: average 1,000 UNIFORM noise readings; repeat 10,000 times
means = rng.uniform(0, 1, (10_000, 1_000)).mean(axis=1)

print(f"μ ≈ {means.mean():.3f}")    # 0.500
print(f"σ ≈ {means.std():.4f}")     # ≈ 0.0091 = 1/√(12·1000)
# histogram `means` and a bell curve appears — from flat ingredients

A single yes/no trial with success probability p — the atom of discrete probability. One coin flip, one transit-or-not observation, one user click. Mean is p, variance is p(1−p) (maximal uncertainty at p=0.5, zero at the extremes). Its direct AI significance: a logistic-regression (or any sigmoid) output IS the p of a Bernoulli — the model says "this light curve is a planet with p=0.83", and binary cross-entropy loss is just the Bernoulli log-likelihood, negated.

• Stack n independent Bernoulli(p) trials and count successes → the Binomial(n, p) distribution.
• The label y ∈ {0, 1} in binary classification is modeled as Bernoulli(p(x)).

import numpy as np
rng = np.random.default_rng(4)

p = 0.3                              # P(transit dip in one observation)
trials = rng.random(10_000) < p      # 10,000 Bernoulli(0.3) draws

print(trials.mean())                 # ≈ 0.30 — mean = p
print(trials.var())                  # ≈ 0.21 — variance = p(1−p)

The arithmetic average — sum divided by count — is the balance point of the data: put the values on a seesaw and the mean is where the fulcrum balances. It's the workhorse summary statistic, but its weakness is the whole reason the next rung exists: a single outlier drags it arbitrarily far, because every value pulls on it with force proportional to its distance.

• The mean minimizes squared distance to all points — which is why MSE-trained models predict means.
• Great for symmetric, well-behaved data; misleading for skewed data (salaries, latencies, planet masses).

import numpy as np

# ten rocky-planet masses (Earth masses)
masses = np.array([0.8, 0.9, 1.0, 1.1, 0.95, 1.05, 0.9, 1.0, 1.1, 0.85])
print(masses.mean())                   # 0.965 — a fair "typical planet"

with_giant = np.append(masses, 317.8)  # one Jupiter sneaks into the sample
print(with_giant.mean())               # 29.77 — "typical" is now nonsense

The middle value when sorted — half the data below, half above. Because it only cares about position in the sorted order, not magnitude, it's robust to outliers where the mean is not: a billionaire joining the room barely moves the median salary but wrecks the mean. When someone says "average", ask which one — for skewed data the choice changes the story completely.

• Median = 50th percentile = p50; the same idea generalizes to p95, p99 — the language of latency SLOs.
• The median minimizes absolute distance to all points (the mean minimizes squared).

import numpy as np

masses = np.array([0.8, 0.9, 1.0, 1.1, 0.95, 1.05, 0.9, 1.0, 1.1, 0.85])
with_giant = np.append(masses, 317.8)   # same billionaire trick, planet edition

print(np.median(masses))       # 0.975
print(np.median(with_giant))   # 1.0   — barely moved
print(with_giant.mean())       # 29.77 — the mean exploded (topic 50)

Variance is the average squared distance from the mean — how spread out the data is. Squaring does two jobs: it kills the sign (deviations below and above both count) and it punishes big deviations disproportionately. The price: variance lives in squared units (km² for km data), which is why its square root — the standard deviation — is what you usually quote.

The n vs n−1 gotcha: dividing by n gives population variance; dividing by n−1 (ddof=1) gives sample variance, which corrects the bias from estimating the mean from the same sample. NumPy defaults to n (ddof=0), pandas defaults to n−1 — a classic source of "why don't our numbers match?" between two scripts.

import numpy as np

periods = np.array([88.0, 224.7, 365.2, 687.0])   # orbital periods (days)

mu = periods.mean()
by_hand = ((periods - mu) ** 2).sum() / len(periods)

print(by_hand == np.var(periods, ddof=0))   # True — population variance (÷n)
print(np.var(periods, ddof=1))              # sample variance (÷(n−1)) — larger

The square root of variance, which puts spread back in the original units — a σ of 0.3 K on a temperature channel means what it says. For Gaussian data, σ comes with the famous 68–95–99.7 rule: about 68% of values fall within 1σ of the mean, 95% within 2σ, 99.7% within 3σ. That rule is the entire basis of "3-sigma anomaly detection": a reading beyond 3σ happens by chance only ~0.3% of the time, so it's probably not chance.

• Z-score: (x − μ) / σ — "how many standard deviations out is this value?"
• Standardizing features = subtract μ, divide by σ — every feature ends up in comparable units.

import numpy as np
rng = np.random.default_rng(5)

noise = rng.normal(0, 1, 100_000)      # telemetry noise, σ = 1

within1 = np.mean(np.abs(noise) < 1)
within2 = np.mean(np.abs(noise) < 2)
within3 = np.mean(np.abs(noise) < 3)
print(f"{within1:.1%} {within2:.1%} {within3:.1%}")   # 68.3% 95.4% 99.7%

The formal way to ask "is this difference real, or just noise?" Frame a null hypothesis H₀ ("no difference"), collect data, and compute how surprising the data would be if H₀ were true. Only if it's very surprising (below a pre-chosen threshold α, usually 0.05) do you reject the null and call the effect significant. This is the machinery behind every A/B test you'll run on models, prompts, and rerankers.

• You never "prove" H₀ — you either reject it or fail to reject it. Absence of evidence ≠ evidence of absence.
• Pick the test for the data: t-test for means, chi-squared for counts, Mann-Whitney when distributions are skewed.
• Decide α and sample size before looking at results — peeking inflates false positives.

import numpy as np
from scipy import stats
rng = np.random.default_rng(6)

old = rng.normal(1.20, 0.3, 200)   # RAG v1 answer latency (s)
new = rng.normal(1.12, 0.3, 200)   # RAG v2 — truly 0.08s faster

t, p = stats.ttest_ind(old, new)
print(f"t={t:.2f}  p={p:.4f}")     # small p → reject "no difference"
# rerun with identical means: p jumps around, mostly large

The p-value is the probability of seeing data at least this extreme if the null hypothesis were true. Small p = "hard to explain as luck." That's all it is — and misreading it is the most common statistical error in tech.

A p-value is NOT the probability the null is true, and 1−p is not the probability your change works. It's a statement about the data under an assumption, not about the hypothesis. And at threshold 0.05, one in twenty truly-null experiments "succeeds" by luck — run enough tests and you will manufacture discoveries. That's the multiple-comparisons problem.

import numpy as np
from scipy import stats
rng = np.random.default_rng(7)

false_alarms = 0
for _ in range(100):                    # 100 A/A tests: no real effect
    a = rng.normal(0, 1, 100)
    b = rng.normal(0, 1, 100)             # identical distributions
    if stats.ttest_ind(a, b).pvalue < 0.05:
        false_alarms += 1
print(false_alarms)                       # ≈ 5 — luck alone "succeeds"

A confidence interval is a range constructed so that, across repeated experiments, it captures the true value 95% of the time (for a 95% CI). It's the honest upgrade from a point estimate: "faithfulness = 0.84" hides everything; "0.84, 95% CI [0.79, 0.89]" tells you the measurement is soft. Note the subtlety: the 95% is a property of the procedure — any single interval either contains the truth or doesn't.

• More data → narrower interval (width shrinks like 1/√n).
• The bootstrap builds a CI with zero formulas: resample your data with replacement, recompute the statistic thousands of times, take the 2.5th and 97.5th percentiles.
• Two models whose CIs overlap heavily? You probably can't tell them apart yet.

import numpy as np
rng = np.random.default_rng(8)

radii = rng.normal(1.6, 0.4, 60)    # 60 measured exoplanet radii (R⊕)

# bootstrap: resample with replacement, recompute the mean 10,000 times
boots = rng.choice(radii, (10_000, 60), replace=True).mean(axis=1)
lo, hi = np.percentile(boots, [2.5, 97.5])
print(f"mean = {radii.mean():.3f}   95% CI = [{lo:.3f}, {hi:.3f}]")

Feature engineering is turning raw columns into signals a model can actually use. Models don't understand your domain — they only see numbers. A raw timestamp is nearly useless; hour, weekday, and is_weekend derived from it can be gold. On tabular problems, one well-crafted feature routinely beats switching to a fancier model, because you're injecting domain knowledge the model could never invent from the raw columns alone.

• Ratios & interactions: mass / radius³ is density — physically meaningful where the raw pair isn't.
• Time components: hour, weekday, month, days-since-event from a single timestamp.
• Aggregations: per-group means, counts, rolling windows ("this star's average planet radius").
• Domain flags: booleans encoding expert rules, like in_habitable_zone.

import pandas as pd
import numpy as np

df = pd.read_csv("exoplanets.csv")  # radius_earths, mass_earths, temp_k, orbital_period

# ratios and transforms that encode physics
df["density"]    = df["mass_earths"] / df["radius_earths"]**3
df["log_period"] = np.log1p(df["orbital_period"])
df["temp_ratio"] = df["temp_k"] / 288          # vs Earth's ~288 K

# a domain flag an astronomer would ask about first
df["earthlike_size"] = df["radius_earths"].between(0.8, 1.6).astype(int)
print(df[["density", "log_period", "temp_ratio", "earthlike_size"]].head())

After engineering comes pruning. Feature selection drops columns that add noise instead of signal — every useless feature is one more dimension to overfit on, one more thing to compute at inference, one more line in the explanation you owe stakeholders. Fewer features means faster, more robust, more explainable models.

Three families, in increasing cost:

• Filter: score each feature alone (correlation with target, mutual information) — cheap, ignores interactions.
• Wrapper: train on subsets and keep the best (RFE) — expensive, respects interactions.
• Embedded: let the model choose while training — tree feature_importances_, L1 coefficients pushed to zero.

from sklearn.ensemble import RandomForestClassifier
import pandas as pd

rf = RandomForestClassifier(n_estimators=300, random_state=0).fit(X_train, y_train)

ranks = pd.Series(rf.feature_importances_, index=X_train.columns)
top5 = ranks.sort_values(ascending=False).head(5).index
print(ranks.sort_values(ascending=False).round(3))

# retrain on just the top 5 and compare
rf_small = RandomForestClassifier(n_estimators=300, random_state=0)
rf_small.fit(X_train[top5], y_train)
print("all:", rf.score(X_test, y_test), "top5:", rf_small.score(X_test[top5], y_test))

Real tables have holes, and most models refuse to train on NaN. The decision is per column, not global: drop the rows (only if few and random), drop the column (if mostly empty), or impute — fill with the median (robust to outliers), the mean, or a per-group value. Fancy option: model the missing value from other columns.

The key insight interviewers fish for: missingness is often a signal itself. A missing mass_earths usually means the planet was found by the transit method, which can't weigh planets — and that correlates with everything else about it. Adding a was_missing 0/1 indicator column lets the model use that fact instead of having it papered over by the median.

from sklearn.impute import SimpleImputer

# add_indicator=True appends the was_missing 0/1 columns automatically
imp = SimpleImputer(strategy="median", add_indicator=True)
X_train_f = imp.fit_transform(X_train)   # fit on train only
X_test_f  = imp.transform(X_test)        # apply the same medians to test

print(X_train.shape, "→", X_train_f.shape)  # extra cols = missingness flags

Data cleaning is everything you do before the data is even wrong-in-interesting-ways: deduplication, type fixes (numbers stored as strings), outlier decisions (is that 40,000 K star a typo or a real O-type?), and normalizing inconsistent categories — 'Kepler', 'KEPLER', and 'kepler ' are three categories to a model and one telescope to you. It's unsexy, and it's where most real-world model failures are born.

The professional habit: cleaning is a script, not a session. Every fix is code that logs what it changed, so the same cleaning runs identically on next month's data pull and an auditor can see exactly what happened to the raw file.

import pandas as pd

df = pd.read_csv("exoplanets_raw.csv")
n0 = len(df)

df = df.drop_duplicates(subset="planet_name")
print(f"dropped {n0 - len(df)} duplicate rows")

# type fix: '3,050' → 3050.0
df["temp_k"] = pd.to_numeric(df["temp_k"].astype(str).str.replace(",", ""), errors="coerce")

# category normalization: 3 spellings → 1
before = df["telescope"].nunique()
df["telescope"] = df["telescope"].str.strip().str.lower()
print(f"telescope categories: {before} → {df['telescope'].nunique()}")

Normalization (min-max scaling) squeezes a feature into [0, 1] with (x − min) / (max − min). The smallest value becomes 0, the largest becomes 1, everything else lands proportionally in between. Simple, interpretable, and exactly what you want when features must share a bounded range — pixel intensities, neural-net inputs, anything fed to a distance metric.

Its Achilles heel: one extreme outlier compresses everyone else into a corner. If orbital periods run 1–400 days and one rogue planet clocks 100,000, then after scaling the rogue sits at 1.0 and every normal planet is crushed into [0, 0.004] — the feature is effectively destroyed. Know your min and max before you trust min-max.

import numpy as np
from sklearn.preprocessing import MinMaxScaler

period = np.array([[3.5], [12.0], [88.0], [225.0], [400.0]])

by_hand = (period - period.min()) / (period.max() - period.min())
sklearn = MinMaxScaler().fit_transform(period)
print(np.allclose(by_hand, sklearn))          # True — same formula

# now inject a rogue 100,000-day orbit and watch the squeeze
rogue = np.vstack([period, [[100_000]]])
print(MinMaxScaler().fit_transform(rogue).round(4).ravel())
# [0. 0.0001 0.0008 0.0022 0.004 1.] — the real planets are gone

Standardization rescales a feature to zero mean and unit variance: z = (x − mean) / std. Unlike min-max it's unbounded, but far more robust to outliers, and it's what distance-based models (KNN, K-Means, SVM) and gradient-based models want — otherwise a feature measured in thousands (temp_k) drowns out one measured in single digits (radius_earths).

The rule that separates professionals from tutorials: fit the scaler on training data only, then apply those same train-derived mean and std to the test set. Fit it on all data and the test set's statistics have leaked into training — your evaluation is quietly optimistic. An sklearn Pipeline enforces this automatically inside cross-validation.

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score

# scaler lives INSIDE the pipeline → refit on each CV training fold only
pipe = make_pipeline(StandardScaler(), KNeighborsClassifier(n_neighbors=5))
scores = cross_val_score(pipe, X, y, cv=5)
print(scores.mean().round(3))

# the leaky version: StandardScaler().fit_transform(X) BEFORE the split —
# test-fold means/stds contaminate training. Never do this.

Models can't do math on strings, so categorical columns must become numbers. The naive fix — kepler=1, tess=2, jwst=3 — invents a fake ordering (is JWST "three times" Kepler?). One-hot encoding avoids that: each category becomes its own 0/1 column, and exactly one is hot per row. No false ordering, and linear models get a separate learnable weight per category.

The breakdown point is high cardinality: 1,000 distinct cities means 1,000 nearly-empty columns — a sparse, memory-hungry matrix where each column carries almost no examples. That's when you reach for target encoding (replace category with its historical target mean) or learned embeddings (Volume 04's territory).

import pandas as pd

df = pd.DataFrame({"telescope": ["kepler", "tess", "jwst", "kepler"],
                   "radius_earths": [1.1, 2.4, 0.9, 11.2]})

encoded = pd.get_dummies(df, columns=["telescope"], dtype=int)
print(encoded.columns.tolist())
# ['radius_earths', 'telescope_jwst', 'telescope_kepler', 'telescope_tess']

# with 10,000 unique values you'd mint 10,000 near-empty columns —
# that's the cue for target encoding or embeddings instead.

Linear regression fits the line (in one dimension) or hyperplane (in many) that minimizes squared error between predictions and targets: y = w·x + b. It's a weak learner but a great teacher: the coefficients are readable ("each extra Earth radius adds 40 K to predicted temperature"), it has a closed-form solution (the normal equation — no iterative training needed), and it's the baseline every fancier model must beat before it earns its complexity.

• Normal equation: w = (XᵀX)⁻¹Xᵀy — solve it directly, no gradient descent.
• Residuals: the vertical gaps between points and the line; squared and summed, they're the loss.
• If a plain line ties your deep model, ship the line.

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
X = rng.uniform(0.2, 5, (200, 1))            # stellar flux (Earth=1)
y = 180 + 95 * X[:, 0] + rng.normal(0, 12, 200)  # temp_k with noise

# normal equation: w = (XᵀX)⁻¹ Xᵀ y  (with a bias column of ones)
Xb = np.c_[np.ones(len(X)), X]
w = np.linalg.solve(Xb.T @ Xb, Xb.T @ y)
print("normal eq :", w.round(2))               # ≈ [180, 95]

lr = LinearRegression().fit(X, y)
print("sklearn   :", lr.intercept_.round(2), lr.coef_.round(2))  # same

Logistic regression is a linear model pushed through a sigmoid, σ(z) = 1/(1+e⁻ᶻ), so the raw score w·x + b becomes a probability between 0 and 1. Despite the name, it's a classifier — "regression" is a historical accident. Above a threshold (default 0.5) you predict positive; below, negative.

It remains the production workhorse for tabular classification baselines: fast to train, calibrated-ish probabilities out of the box, coefficients that show which features push toward which class, and hard to overfit with a little regularization. The probability output matters as much as the label — 0.93 and 0.51 are both "positive," but you'd only wake someone up for the first.

from sklearn.linear_model import LogisticRegression

# X: [temp_ratio, radius_earths, density] · y: 1 = habitable-zone candidate
clf = LogisticRegression(max_iter=1000).fit(X_train, y_train)

proba = clf.predict_proba(X_test)[:, 1]   # P(candidate) per planet
print(proba[:4].round(3))                 # e.g. [0.93 0.08 0.51 0.77]
print(clf.predict(X_test)[:4])            # same rows thresholded at 0.5

# which features push toward 'candidate'?
print(dict(zip(X_train.columns, clf.coef_[0].round(2))))

A decision tree is recursive if/else: at each node, pick the feature and threshold whose split makes the two child groups purest (measured by Gini impurity or entropy gain), then recurse into each child. Prediction is just walking a path from root to leaf — which means the model is a human-readable ruleset.

The catch: grown to full depth, a tree will happily memorize the training set — one leaf per weird example, perfect train accuracy, poor generalization. Depth limits and minimum-leaf sizes tame it, but the real move is the next three items: trees are weak alone but powerful in ensembles, which is why the forest and boosting rungs come right after.

from sklearn.tree import DecisionTreeClassifier, export_text

tree = DecisionTreeClassifier(max_depth=3, random_state=0).fit(X_train, y_train)
print(export_text(tree, feature_names=list(X_train.columns)))
# |--- temp_k <= 350.0
# |   |--- radius_earths <= 1.8
# |   |   |--- class: candidate ...

# unlimited depth memorizes: near-perfect train, worse test
deep = DecisionTreeClassifier(random_state=0).fit(X_train, y_train)
print(deep.score(X_train, y_train), deep.score(X_test, y_test))

A random forest fixes the single tree's memorization problem with deliberate diversity: train many trees, each on a bootstrap sample of the rows (drawn with replacement) and considering only a random subset of features at each split. Then average their votes. Each tree still overfits — but each overfits differently, and the errors largely cancel when you average. That's variance reduction by ensembling (this technique is called bagging).

• The reliable default for tabular data: strong out of the box with almost no tuning.
• Parallel: trees are independent, so training scales across cores.
• Free extras: feature_importances_ and out-of-bag error estimates.

from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier

single = DecisionTreeClassifier(random_state=0).fit(X_train, y_train)
forest = RandomForestClassifier(n_estimators=200, random_state=0).fit(X_train, y_train)

print("single tree :", single.score(X_test, y_test).round(3))
print("200-tree RF :", forest.score(X_test, y_test).round(3))
# same data, same split — the forest usually wins by a clear margin

Where the forest builds trees in parallel and averages, gradient boosting builds them sequentially: each new tree is trained to fix the errors (residuals) of the ensemble so far, and its prediction is added with a small learning rate. Dozens to thousands of small corrections compound into a very strong model. XGBoost is the implementation that ruled tabular competitions for a decade: strong regularization built into the objective, second-order gradients, sparsity-aware splits, fast histogram training.

The knob that matters most in practice: early stopping. Watch validation loss as rounds accumulate and stop when it stops improving — the model picks its own size instead of you guessing n_estimators.

from xgboost import XGBClassifier

clf = XGBClassifier(
    n_estimators=2000,           # generous cap — early stopping decides
    learning_rate=0.05,
    max_depth=5,
    early_stopping_rounds=50,    # stop after 50 rounds with no val gain
    eval_metric="auc",
)
clf.fit(X_train, y_train, eval_set=[(X_val, y_val)], verbose=False)
print("stopped at round:", clf.best_iteration)
print("val AUC:", clf.best_score)

LightGBM is gradient boosting rebuilt for speed on large data. Two ideas do the heavy lifting. First, histogram binning: continuous features are pre-bucketed into ~255 bins, so finding the best split means scanning bins, not sorting raw values. Second, leaf-wise growth: instead of growing every tree level by level (as XGBoost classically did), LightGBM always splits the leaf with the biggest loss reduction — deeper, more asymmetric trees that get more accuracy per split.

The trade: leaf-wise trees can overfit small datasets, so num_leaves and min_data_in_leaf are the knobs to respect. On big tabular data, it's usually the fastest strong option — often several times quicker than XGBoost at equal accuracy.

import time
from lightgbm import LGBMClassifier
from xgboost import XGBClassifier
from sklearn.metrics import roc_auc_score

for name, model in [("lgbm", LGBMClassifier(n_estimators=500)),
                    ("xgb ", XGBClassifier(n_estimators=500))]:
    t0 = time.perf_counter()
    model.fit(X_train, y_train)
    auc = roc_auc_score(y_test, model.predict_proba(X_test)[:, 1])
    print(f"{name}: {time.perf_counter()-t0:5.1f}s  AUC={auc:.4f}")

CatBoost is the boosting library that eats categorical features natively — you pass raw strings and name the columns; no one-hot step. Internally it uses ordered target statistics: each row's category gets encoded using the target mean of only the rows that came before it in a random ordering. That "only the past" trick prevents the row's own label from leaking into its own encoding — the classic failure of naive target encoding, and the reason CatBoost overfits less on high-cardinality columns.

• Great defaults: often competitive with tuned XGBoost/LightGBM untouched.
• Also uses ordered boosting and symmetric trees for extra overfitting resistance.
• The pick when your table is full of string columns like facilities, spectral types, instrument names.

from catboost import CatBoostClassifier

# raw string columns go straight in — no get_dummies, no target-encode step
cat_cols = ["telescope", "discovery_facility", "spectral_type"]

clf = CatBoostClassifier(iterations=800, learning_rate=0.05, verbose=200)
clf.fit(X_train, y_train,
        cat_features=cat_cols,
        eval_set=(X_val, y_val))
print(clf.get_best_score())

A support vector machine asks a sharper question than "find a separating boundary": find the boundary with maximum margin — the greatest distance to the nearest training points of either class. Those nearest points are the support vectors; they alone define the boundary, and every other point could be deleted without changing it. Fat margins generalize better than boundaries that graze the data.

The second idea is the kernel trick: to get a non-linear boundary, implicitly map points to a higher-dimensional space where they are linearly separable — without ever computing that mapping. The RBF kernel does this with a similarity function, turning "line" into "any smooth curved boundary." Scale your features first (item 62); SVMs are distance-based to the core.

from sklearn.datasets import make_moons
from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

X, y = make_moons(n_samples=400, noise=0.25, random_state=0)

for kernel in ("linear", "rbf"):
    pipe = make_pipeline(StandardScaler(), SVC(kernel=kernel))
    print(kernel, pipe.fit(X, y).score(X, y).round(3))
# linear ≈ 0.88 (a line can't bend around moons) · rbf ≈ 0.97

K-nearest neighbors has no training step at all: to classify a point, find the k closest training points (usually Euclidean distance) and take a majority vote. The "model" is literally the training data. That makes it trivially simple, but expensive at inference (every prediction scans the dataset) and cursed by dimensionality — in high dimensions, everything is roughly equidistant from everything, and "nearest" loses meaning.

Why it matters beyond the basics: KNN is the conceptual ancestor of vector search. Retrieving the top-k most similar embeddings from a vector database is exactly a nearest-neighbor query — the same idea, industrialized with approximate indexes (HNSW, IVF) so it scales to billions of vectors.

import numpy as np
from sklearn.neighbors import KNeighborsClassifier

def knn3(X_train, y_train, q):
    dists = np.linalg.norm(X_train - q, axis=1)   # Euclidean to every row
    top3 = np.argsort(dists)[:3]                   # indices of 3 nearest
    return np.bincount(y_train[top3]).argmax()     # majority vote

mine = [knn3(X_tr, y_tr, q) for q in X_te]
skl  = KNeighborsClassifier(n_neighbors=3).fit(X_tr, y_tr).predict(X_te)
print((mine == skl).all())   # True — same votes

No labels? K-Means finds structure anyway. Pick k centroids, then loop two steps until nothing moves: assign every point to its nearest centroid, then move each centroid to the mean of its assigned points. It converges fast and runs everywhere — but you must choose k yourself, and it assumes clusters are roughly round and similar in size (it partitions space into nearest-centroid cells, which can't wrap around crescents or chains).

• Inertia: sum of squared point-to-centroid distances — K-Means' internal loss.
• Elbow method: plot inertia for k = 1…10; pick k where the curve stops dropping steeply.
• Standardize features first — it's pure distance math.

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

Xs = StandardScaler().fit_transform(planets[["log_period", "radius_earths", "temp_k"]])

for k in range(1, 11):
    km = KMeans(n_clusters=k, n_init="auto", random_state=0).fit(Xs)
    print(k, round(km.inertia_, 1))   # look for the elbow

labels = KMeans(n_clusters=4, n_init="auto", random_state=0).fit_predict(Xs)

DBSCAN clusters by density, not distance-to-centroid: a point with at least min_samples neighbors within radius eps is a core point; cores within reach of each other chain into clusters of any shape — crescents, rings, filaments. Points in sparse regions join nothing and are labeled noise (label −1). Two superpowers K-Means lacks: no k to choose, and an explicit outlier concept built in.

The cost: you tune eps and min_samples instead, and clusters with very different densities are hard to capture with a single eps. Think of DBSCAN as the complement to K-Means — reach for it when shapes are weird or when the noise points are what you actually care about.

from sklearn.cluster import DBSCAN, KMeans
from sklearn.datasets import make_moons
import numpy as np

X, _ = make_moons(n_samples=400, noise=0.07, random_state=0)

db = DBSCAN(eps=0.2, min_samples=5).fit(X)
print("clusters:", len(set(db.labels_)) - (1 if -1 in db.labels_ else 0))
print("noise pts:", np.sum(db.labels_ == -1))
# DBSCAN: 2 clean moons · KMeans(2) slices straight across both

Hierarchical (agglomerative) clustering doesn't pick a number of clusters — it builds every clustering at once. Start with each point as its own cluster; repeatedly merge the two closest clusters until one remains. The merge history forms a dendrogram: a tree whose branch heights record how dissimilar the merged groups were. Cut the tree at any height and you get a clustering; lower cuts give many small clusters, higher cuts give few large ones.

Often the dendrogram itself is the deliverable — it shows structure at every scale simultaneously, which no single flat clustering can. The main knob is linkage: how to measure distance between clusters (ward minimizes variance and is the usual default; single chains; complete compacts).

from scipy.cluster.hierarchy import linkage, dendrogram, fcluster
import matplotlib.pyplot as plt

Z = linkage(Xs, method="ward")      # the full merge tree
dendrogram(Z, truncate_mode="level", p=5)
plt.show()

# cut at two different heights — no refit needed
coarse = fcluster(Z, t=2, criterion="maxclust")   # 2 clusters
fine   = fcluster(Z, t=6, criterion="maxclust")   # 6 clusters, same tree

Principal Component Analysis finds the directions along which your data varies the most — the eigenvectors of the covariance matrix (rung 36 returns) — and projects the data onto the top few. PC1 is the single axis capturing maximum variance; PC2 is the best axis perpendicular to it; and so on. Keep 2 components of a 64-dimensional dataset and you have the best possible 2D "shadow" of it, with explained_variance_ratio_ telling you exactly how much information the shadow retains.

• Compression: 64 features → 10 components that keep 95% of the variance.
• First look: the standard way to eyeball high-dimensional data in 2D.
• Linear only: it can only rotate and project — curved structure needs t-SNE/UMAP (next two items).

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

Xs = StandardScaler().fit_transform(X)      # PCA needs comparable scales
pca = PCA(n_components=2)
X2 = pca.fit_transform(Xs)

print(pca.explained_variance_ratio_)        # e.g. [0.46 0.21] → 67% kept
plt.scatter(X2[:, 0], X2[:, 1], c=y, s=8, cmap="cool")
plt.xlabel("PC1"); plt.ylabel("PC2"); plt.show()

t-SNE is a non-linear projection built for one job: revealing cluster structure. It converts distances into neighbor probabilities in high dimensions, then arranges points in 2D so that local neighborhoods are preserved — points that were close stay close. Where PCA smears curved manifolds flat, t-SNE untangles them into crisp islands.

Now the warning label, which interviewers love: in a t-SNE plot, distances between clusters mean nothing and cluster sizes mean nothing. The algorithm expands dense regions and compresses sparse ones, and different perplexity values redraw the map. Read it as "these points group together" — never as "cluster A is twice as far from B as from C."

from sklearn.datasets import load_digits
from sklearn.manifold import TSNE
import matplotlib.pyplot as plt

digits = load_digits()                      # 64-dimensional images
X2 = TSNE(n_components=2, perplexity=30,
          random_state=0).fit_transform(digits.data)

plt.scatter(X2[:, 0], X2[:, 1], c=digits.target, s=6, cmap="tab10")
plt.title("10 digit islands — compare with your PCA plot"); plt.show()

UMAP plays the same game as t-SNE — non-linear projection preserving local neighborhoods — but built on manifold theory, it's typically much faster, scales to hundreds of thousands of points, and preserves global structure better: the relative placement of islands is somewhat more meaningful (though still not literal distances). It also supports .transform() on new points, which t-SNE classically can't.

For an AI engineer, UMAP's killer app is eyeballing embedding spaces: project your vector database's embeddings to 2D and you can literally see whether topics form clean islands, where duplicate chunks pile up, and which queries land in no-man's-land. Key knobs: n_neighbors (small = local detail, large = global shape) and min_dist (how tightly points pack).

import umap
from sentence_transformers import SentenceTransformer
import matplotlib.pyplot as plt

model = SentenceTransformer("all-MiniLM-L6-v2")
emb = model.encode(headlines)               # 200 space-news headlines → 384-dim

X2 = umap.UMAP(n_neighbors=15, min_dist=0.1,
               random_state=0).fit_transform(emb)
plt.scatter(X2[:, 0], X2[:, 1], s=10)
plt.title("do topics form islands?"); plt.show()

Accuracy is the fraction of predictions that were correct: (TP + TN) / total. It's the first metric everyone learns and the one most likely to lie to you. On imbalanced data it's deceptive to the point of uselessness: if 99% of rows are negative, a model that always predicts negative scores 99% accuracy while catching exactly zero positives — the entire reason the next five rungs of this ladder exist.

• Trust it only when classes are roughly balanced and all errors cost the same.
• Always compare against the majority-class baseline: if "always predict the biggest class" scores 0.99, your 0.99 model has learned nothing.

from sklearn.dummy import DummyClassifier
from sklearn.metrics import accuracy_score, recall_score

# y_test: 990 zeros, 10 ones (habitable-zone candidates)
lazy = DummyClassifier(strategy="most_frequent").fit(X_train, y_train)
pred = lazy.predict(X_test)

print("accuracy:", accuracy_score(y_test, pred))   # 0.99 — looks amazing
print("recall  :", recall_score(y_test, pred))     # 0.00 — caught nothing

Precision answers: of everything the model flagged positive, how much really was? precision = TP / (TP + FP). It's the metric of false-alarm cost. Optimize precision when acting on a flag is expensive or damaging: a spam filter that junks real mail, a fraud system that blocks legitimate cards, an alert queue that pages a human.

Precision looks only at the flagged pile — it says nothing about the positives you never flagged (that's recall, next item). You can make precision arbitrarily high by flagging only your single most confident case, which is exactly why precision and recall are always read together.

from sklearn.metrics import precision_score, confusion_matrix

# imaginary spam filter: write out the matrix, compute by hand first
#                 predicted spam   predicted ham
# actual spam          TP=40           FN=10
# actual ham           FP=5            TN=945
tp, fp = 40, 5
print("by hand :", tp / (tp + fp))                 # 0.889

print("sklearn :", precision_score(y_test, model.predict(X_test)).round(3))

Recall answers the mirror question: of everything actually positive, how much did the model catch? recall = TP / (TP + FN). It's the metric of miss cost. Optimize recall when a missed positive is the disaster: cancer screening, safety-content filters, structural-fault detection. A screening test would rather raise ten false alarms than let one real case walk out the door — the alarms get a second look; the miss never does.

Precision and recall trade off through the decision threshold: lower it and you flag more (recall up, precision down); raise it and you flag less (precision up, recall down). Which side you favor is a product decision about which error hurts more — not a modeling decision.

from sklearn.metrics import recall_score

# same spam-filter matrix as item 80
tp, fn = 40, 10
print("by hand :", tp / (tp + fn))                 # 0.800

# lower the threshold → flag more → recall rises, precision falls
proba = model.predict_proba(X_test)[:, 1]
for t in (0.8, 0.5, 0.2):
    print(t, recall_score(y_test, proba >= t).round(3))

When you need one number that respects both precision and recall, the F1 score is their harmonic mean: F1 = 2PR / (P + R). Harmonic — not arithmetic — and that choice has teeth: the harmonic mean is dragged toward the smaller value, so a model can't buy a good F1 by maxing one metric and abandoning the other. Precision 0.9 with recall 0.3 averages to 0.6 arithmetically, but F1 says 0.45 — it punishes the imbalance hard.

• F1 = 1 only when precision and recall are both 1; either at 0 makes F1 exactly 0.
• Use it to compare models on imbalanced data with one scalar; use the P/R pair when the two error costs differ a lot.
• Variants: F0.5 weights precision more, F2 weights recall more.

from sklearn.metrics import f1_score

p, r = 0.9, 0.3
f1 = 2 * p * r / (p + r)
print("arithmetic mean:", (p + r) / 2)   # 0.60 — flattering
print("F1 (harmonic)  :", round(f1, 3))  # 0.45 — honest

# on real predictions
print(f1_score(y_test, model.predict(X_test)).round(3))

Every metric so far depends on one chosen threshold. ROC-AUC judges the model across all thresholds at once: sweep the cutoff from 1 down to 0, and at each point plot true-positive rate (recall) against false-positive rate. That trace is the ROC curve; the area under it summarizes ranking quality. AUC has a beautiful interpretation: it's the probability that a randomly chosen positive gets a higher score than a randomly chosen negative.

• 0.5 = coin-flip ranking (the diagonal); 1.0 = every positive ranked above every negative.
• Threshold-free: perfect for comparing models before you commit to an operating point.
• On extreme imbalance, precision-recall AUC is often more informative — FPR barely moves when negatives dominate.

from sklearn.metrics import roc_curve, roc_auc_score
import matplotlib.pyplot as plt

for name, m in [("forest", rf), ("xgboost", xgb)]:
    proba = m.predict_proba(X_test)[:, 1]
    fpr, tpr, _ = roc_curve(y_test, proba)
    auc = roc_auc_score(y_test, proba)
    plt.plot(fpr, tpr, label=f"{name} AUC={auc:.3f}")

plt.plot([0, 1], [0, 1], "--", label="chance")
plt.legend(); plt.show()

The confusion matrix is the 2×2 grid every metric above is computed from: rows are the truth, columns are the prediction, and the four cells are TP (hit), FP (false alarm), FN (miss), and TN (correct rejection). Accuracy is the diagonal over the total; precision is the left column's purity; recall is the top row's coverage. Reading one fluently — including instantly spotting which error type dominates — is table stakes in any ML conversation.

The habit to build: look at the matrix before any scalar metric. Two models with identical F1 can fail in opposite directions — one drowning in false alarms, the other quietly missing positives — and only the matrix shows which. For multi-class problems it generalizes to k×k, where off-diagonal hot spots reveal exactly which classes get confused with which.

from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay
import matplotlib.pyplot as plt

cm = confusion_matrix(y_test, model.predict(X_test))
tn, fp, fn, tp = cm.ravel()      # sklearn order: [[TN, FP], [FN, TP]]
print(f"hits {tp} · misses {fn} · false alarms {fp} · correct rejects {tn}")

ConfusionMatrixDisplay(cm, display_labels=["non-hab", "candidate"]).plot(cmap="Blues")
plt.show()

PyTorch is built from three ideas that stack cleanly. A tensor is a NumPy array that can live on a GPU. Autograd watches every operation you perform on tensors that have requires_grad=True and silently records a computation graph — so when you call loss.backward(), it walks that graph in reverse and fills in .grad for every parameter automatically. nn.Module packages weights and a forward() method into composable layers. That's the whole framework; everything else is convenience.

The reason PyTorch won research and then production is that the training loop is your Python code — an ordinary for loop you can print from, debug, and modify — not a black box. Every LLM you'll touch later (Volume 05 onward) is a big nn.Module trained by exactly the loop below.

• The sacred loop: forward → loss → backward() → step() → zero_grad(). Memorize the order.
• Forgetting zero_grad() makes gradients accumulate across steps — the classic first bug.
• .to("cuda") moves tensors/models to GPU; model and data must be on the same device.

import torch

# learn y = 2x + 1 from scratch: two parameters, the sacred loop
x = torch.linspace(0, 10, 50).unsqueeze(1)
y = 2 * x + 1

w = torch.randn(1, requires_grad=True)
b = torch.zeros(1, requires_grad=True)
opt = torch.optim.SGD([w, b], lr=0.01)

for epoch in range(200):
    loss = ((w * x + b - y) ** 2).mean()  # forward + MSE
    loss.backward()                          # autograd fills w.grad, b.grad
    opt.step()                               # w -= lr * w.grad (roughly)
    opt.zero_grad()

print(w.item(), b.item())  # → ≈2.0, ≈1.0

TensorFlow is Google's framework, and its personality is the opposite of PyTorch's: it wants to compile your model into a graph. Modern TF runs eagerly by default (line-by-line, like PyTorch), but decorate a function with @tf.function and TF traces it once, freezes it into an optimized dataflow graph, and reruns that graph without touching Python again — great for speed, awkward for debugging (your print() fires only during tracing).

You'll likely write PyTorch, but you must read TensorFlow fluently: an enormous amount of production ML — especially older pipelines, mobile apps (TF Lite), and serving infrastructure (TF Serving, the SavedModel format) — is TF. Interviewers also love the "eager vs graph" question because it tests whether you understand what a framework actually does.

• Eager mode: ops execute immediately — TF 2.x default, feels like NumPy.
• @tf.function: traces Python into a graph; Python side effects happen once, tensor ops every call.
• SavedModel: a directory format bundling graph + weights — the deployment contract for TF Serving.

import tensorflow as tf

w = tf.Variable(0.0); b = tf.Variable(0.0)
opt = tf.keras.optimizers.SGD(learning_rate=0.01)

@tf.function   # ← everything below is traced into a graph on first call
def train_step(x, y):
    with tf.GradientTape() as tape:      # TF's autograd is explicit
        loss = tf.reduce_mean((w * x + b - y) ** 2)
    grads = tape.gradient(loss, [w, b])
    opt.apply_gradients(zip(grads, [w, b]))
    return loss

x = tf.linspace(0.0, 10.0, 50); y = 2 * x + 1
for _ in range(200): train_step(x, y)
print(w.numpy(), b.numpy())   # → ≈2.0, ≈1.0

Keras is the high-level API that turns "define, compile, fit" into three declarative calls. You describe the architecture as a stack of layers, model.compile() attaches an optimizer, loss, and metrics, and model.fit() runs the entire training loop — batching, shuffling, validation splits, progress bars — that you wrote by hand in #85. Since Keras 3 it's multi-backend: the same code can run on TensorFlow, JAX, or PyTorch underneath.

The tradeoff is honest: you gain speed for standard architectures and lose sight of the loop. Use Keras when the problem is conventional and the value is iteration speed; drop to raw PyTorch when you need a custom loss, a weird sampling scheme, or gradient surgery. Knowing what fit() hides — because you built it yourself in #85 — is what separates an engineer from an API user.

import keras
from keras import layers

# 3-layer classifier: 200 flux features → transit / no-transit
model = keras.Sequential([
    layers.Dense(64, activation="relu", input_shape=(200,)),
    layers.Dense(32, activation="relu"),
    layers.Dense(1,  activation="sigmoid"),
])
model.compile(optimizer="adam", loss="binary_crossentropy",
              metrics=["accuracy"])
model.fit(X_train, y_train, epochs=10, batch_size=32,
          validation_split=0.2)          # the whole #85 loop, one line

The perceptron is one artificial neuron: multiply each input by a learned weight, add them up with a bias, and pass the sum through a hard threshold — output 1 if the sum is positive, else 0. Geometrically it draws a single straight line (a hyperplane) through input space and classifies by which side you're on. Its learning rule is beautifully simple: when it misclassifies, nudge the weights toward the correct answer; Rosenblatt proved this converges whenever a separating line exists.

Its fame comes from its failure: XOR — output 1 when inputs differ — has no separating line, so a single perceptron can literally never learn it. Minsky and Papert's 1969 proof of this froze neural-net research for a decade, and the fix (stack layers, #89) is the reason "deep" learning exists. Every interview question about "why do we need hidden layers?" is really asking about this.

import numpy as np

X = np.array([[0,0],[0,1],[1,0],[1,1]])
y = np.array([0, 0, 0, 1])          # AND — linearly separable
w, b, lr = np.zeros(2), 0.0, 0.1

for epoch in range(10):
    for xi, yi in zip(X, y):
        pred = int(xi @ w + b > 0)
        w += lr * (yi - pred) * xi     # the perceptron learning rule
        b += lr * (yi - pred)

print([int(xi @ w + b > 0) for xi in X])  # [0,0,0,1] ✓
# swap y for XOR = [0,1,1,0] and this loop cycles forever — no line exists

A neural network is layers of perceptron-like units where the hard threshold is replaced by a smooth non-linear activation (ReLU, tanh, sigmoid), stacked so each layer's output feeds the next. The non-linearity is not decoration — it is the entire point. Compose two linear maps and you get… another linear map: W₂(W₁x) = (W₂W₁)x. A hundred linear layers collapse into one matrix. Insert a non-linearity between them and the collapse is broken: the network can now bend, fold, and warp input space.

With enough hidden units, one hidden layer can approximate any continuous function — the universal approximation theorem. That's why the same recipe (linear → activation → linear → …) classifies light curves, translates languages, and powers every LLM: it's a general-purpose function-fitting machine, and the rest of this volume is about how to fit it efficiently.

• Width = units per layer (capacity per step); depth = number of layers (composition of steps).
• ReLU max(0, x) is the modern default: cheap, and its gradient doesn't vanish for positive inputs.
• Remove the activations and any network — however deep — degrades to linear regression.

import torch, torch.nn as nn

X = torch.tensor([[0.,0.],[0.,1.],[1.,0.],[1.,1.]])
y = torch.tensor([[0.],[1.],[1.],[0.]])       # XOR — the perceptron-killer

model = nn.Sequential(nn.Linear(2, 3), nn.Tanh(),   # delete Tanh → never learns
                      nn.Linear(3, 1), nn.Sigmoid())
opt = torch.optim.Adam(model.parameters(), lr=0.1)

for _ in range(500):
    loss = nn.functional.binary_cross_entropy(model(X), y)
    opt.zero_grad(); loss.backward(); opt.step()

print(model(X).detach().round().squeeze())  # → [0., 1., 1., 0.] ✓

Forward propagation is the network actually computing: data enters, and each layer performs one matrix multiply, one bias add, one activation — a = activation(W·x + b) — then hands the result to the next layer. That's it. Inference is one forward pass; training is a forward pass plus a backward pass (#91). Every intimidating architecture diagram you'll ever see — ResNet, Transformer, diffusion U-Net — is just a fancier routing of this same per-layer computation.

Hand-computing one full forward pass, once, with real numbers, is the highest-value exercise in this volume. It makes tensor shapes concrete (a (2,3) weight matrix turns a 2-vector into a 3-vector), and shape errors are the single most common deep-learning bug you'll debug.

import numpy as np

# fixed weights → verify your paper computation to 4 decimals
W1 = np.array([[0.5, -0.3], [1.0, 0.8]]); b1 = np.array([0.0, 0.0])
W2 = np.array([[1.2, -0.7]]);              b2 = np.array([0.1])

x  = np.array([1.0, 0.0])
z1 = W1 @ x + b1          # [0.5, 1.0]   — linear
h  = np.tanh(z1)          # [0.4621, 0.7616] — bend
z2 = W2 @ h + b2          # [0.1214]
y  = 1 / (1 + np.exp(-z2))
print(h.round(4), y.round(4))   # [0.4621 0.7616] [0.6303]

This is the chain rule from calculus (topic #40) industrialized: if L depends on ŷ, which depends on h, which depends on W₁, then ∂L/∂W₁ = ∂L/∂ŷ · ∂ŷ/∂h · ∂h/∂W₁. Backprop is nothing more — and nothing less — than organizing that multiplication so no term is ever computed twice.

Backpropagation answers the only question training ever asks: "for each of the network's million weights, if I nudged it slightly, how much would the loss change?" That number is the weight's gradient, and the naive way to get it — nudge each weight, rerun the forward pass, compare — would cost one million forward passes per step. Backprop gets every gradient in a single backward sweep, for roughly the cost of one extra forward pass. That efficiency is the entire reason deep learning is computationally possible.

The mechanism: during the forward pass, each layer caches its inputs. Then, starting from the loss, an error signal ∂L/∂output flows backward. At each layer the signal does two jobs: (1) multiplied by the cached input, it yields that layer's weight gradients — which get stored for the optimizer; (2) multiplied by the layer's weights, it becomes the error signal for the layer below — and flows on. Every layer only needs local information: its cached input, its own weights, and the signal arriving from above. No layer knows anything about the whole network — yet the chain rule guarantees the global gradients come out exactly right.

This is the crux rung of all 308 topics. Vanishing/exploding gradients (#98), why BatchNorm helps (#96), why ResNet skip connections exist (#100), what Adam consumes (#94), what LoRA computes gradients through, why long sequences blow up transformer memory — every one of these is a statement about backprop's error signal. Own this and the rest of the ladder gets easier.

• Forward caches, backward consumes: that's why training uses far more memory than inference.
• Each layer multiplies the incoming signal by its local derivative — a chain of factors. Many factors < 1 → signal vanishes; many > 1 → it explodes.
• loss.backward() in PyTorch runs precisely this sweep over the recorded autograd graph.
• Cost: one backward sweep ≈ 2× a forward pass, regardless of parameter count.

import numpy as np; import torch

# tiny net: y_hat = sigmoid(w2 · tanh(w1 · x)) — backprop by hand
w1, w2, x, y = 0.5, -0.4, 1.5, 1.0

h     = np.tanh(w1 * x)                    # forward, cache everything
y_hat = 1 / (1 + np.exp(-w2 * h))
L     = (y_hat - y) ** 2

dL_dyhat = 2 * (y_hat - y)                 # backward: chain rule, right to left
dyhat_dz = y_hat * (1 - y_hat)             # sigmoid's local derivative
dL_dw2   = dL_dyhat * dyhat_dz * h         # drop gradient for w2
dL_dh    = dL_dyhat * dyhat_dz * w2        # signal keeps flowing down…
dL_dw1   = dL_dh * (1 - h ** 2) * x        # tanh' = 1 − h², drop w1's gradient

# verify against autograd — the numbers must match exactly
tw1 = torch.tensor(w1, requires_grad=True); tw2 = torch.tensor(w2, requires_grad=True)
tl = (torch.sigmoid(tw2 * torch.tanh(tw1 * x)) - y) ** 2
tl.backward()
print(dL_dw1, tw1.grad.item())   # same number, twice ✓
print(dL_dw2, tw2.grad.item())   # same number, twice ✓

Backprop (#91) gives you every weight's gradient — the direction of steepest loss increase. Gradient descent is the obvious response: step every weight the opposite way, w ← w − lr·∂L/∂w, and repeat until the loss stops falling. The learning rate lr scales the step, and it is the most consequential hyperparameter in all of deep learning.

Picture the loss as a valley. Too large a learning rate and you leap across the valley, land higher than you started, leap again — the loss diverges to NaN. Too small and you inch downhill for hours, possibly stalling on a flat shelf. The classic tell: loss shooting up or oscillating → cut lr by 10×; loss creeping down a fraction of a percent per epoch → raise it. Modern training wraps lr in schedules (warmup, cosine decay), but the intuition here is what those schedules are managing.

import torch

def train(lr, steps=60):
    w = torch.tensor([8.0], requires_grad=True)
    for _ in range(steps):
        loss = (w ** 2).sum()          # a simple bowl, minimum at w=0
        loss.backward()
        with torch.no_grad(): w -= lr * w.grad
        w.grad.zero_()
    return loss.item()

for lr in [10, 0.1, 0.0001]:
    print(f"lr={lr}: final loss {train(lr):.2e}")
# lr=10 → inf (diverged) · lr=0.1 → ~0 (converged) · lr=1e-4 → 62.9 (barely moved)

True gradient descent computes the gradient over the entire dataset before taking one step — perfectly accurate, absurdly expensive at scale. SGD instead estimates the gradient from a small random minibatch (32, 64, 256 samples), takes a slightly-wrong step, and repeats. You take thousands of cheap noisy steps instead of a few expensive perfect ones, and the total progress per unit of compute is dramatically better.

The surprise is that the noise isn't just tolerable — it helps. Noisy steps jitter the trajectory out of sharp, narrow minima (which generalize poorly) and settle it into wide, flat basins (which generalize well). This is why minibatch SGD-family optimizers train literally every modern model, and why batch size is a real hyperparameter: bigger batches = smoother but more expensive steps and less of that useful regularizing noise.

• Epoch = one full pass over the data; with batch 32 and 90k samples that's ~2,800 steps per epoch.
• Shuffle every epoch — otherwise the "random" batches repeat and the noise gains a pattern.
• In PyTorch, "SGD" = torch.optim.SGD fed by a DataLoader; the loader provides the stochasticity.

import torch
from torch.utils.data import TensorDataset, DataLoader

ds = TensorDataset(flux, labels)          # 90k light curves
loader = DataLoader(ds, batch_size=32, shuffle=True)  # ← the "stochastic"
opt = torch.optim.SGD(model.parameters(), lr=0.01, momentum=0.9)

for epoch in range(5):
    for xb, yb in loader:                 # ~2,800 noisy steps per epoch
        loss = loss_fn(model(xb), yb)     # gradient from 32 samples only
        opt.zero_grad(); loss.backward(); opt.step()
# full-batch equivalent: ONE step per epoch, each needing all 90k forward passes

Plain SGD gives every one of a million weights the same learning rate. But real networks are lopsided: some weights see huge, steady gradients; others see rare, tiny ones. Adam keeps two exponential moving averages per parameter — m, the mean of recent gradients (momentum: where have gradients been pointing?), and v, the mean of recent squared gradients (scale: how big have they been?) — and steps each weight by lr · m̂ / (√v̂ + ε). Consistently large gradients get damped; small steady ones get amplified. Every parameter effectively tunes its own learning rate.

In practice you'll use AdamW, which fixes a subtle flaw: Adam entangles L2 weight decay with its adaptive scaling, so decay strength varies per parameter. AdamW decouples the decay — applying it directly to weights, outside the adaptive machinery — and it is the default optimizer for pretraining and fine-tuning essentially every transformer and LLM you will ever touch. Default starting point: lr=3e-4, or 1e-5–5e-5 when fine-tuning pretrained models.

# the only line that changes vs SGD — same loop, same everything else
opt = torch.optim.AdamW(model.parameters(), lr=3e-4, weight_decay=0.01)

for epoch in range(5):
    for xb, yb in loader:
        loss = loss_fn(model(xb), yb)
        opt.zero_grad(); loss.backward(); opt.step()

# typical outcome on the same net: SGD lr=0.01 needs ~40 epochs to reach
# the loss AdamW hits in ~8 — per-parameter adaptation does the tuning for you

RMSProp is Adam's direct ancestor and the cleanest way to understand what "adaptive learning rate" means. It keeps just one running average per parameter — the RMS (root-mean-square) of recent gradients — and divides the learning rate by it: w ← w − lr · g / (RMS(g) + ε). Parameters whose gradients have been large get small careful steps; parameters with faint gradients get boosted ones. Geoff Hinton introduced it in a Coursera lecture (it was never formally published), and it fixed AdaGrad's fatal flaw of accumulating squared gradients forever until learning ground to a halt.

Hold the family tree in your head and optimizer questions become easy: SGD is the baseline; add a velocity term and you get momentum; divide by the gradient RMS and you get RMSProp; do both (plus bias correction for the early steps) and you get Adam (#94). You'll still meet RMSProp in the wild — classic reinforcement-learning configs (DQN) and older Keras defaults — and knowing it as "Adam minus momentum" means those codebases read instantly.

import torch

optimizers = {
    "sgd":     lambda p: torch.optim.SGD(p, lr=0.01),
    "rmsprop": lambda p: torch.optim.RMSprop(p, lr=0.01, alpha=0.99),
    "adam":    lambda p: torch.optim.Adam(p, lr=0.01),
}
for name, make in optimizers.items():
    model = build_model(); opt = make(model.parameters())
    losses = []
    for step, (xb, yb) in enumerate(loader):
        loss = loss_fn(model(xb), yb)
        opt.zero_grad(); loss.backward(); opt.step()
        losses.append(loss.item())
    print(name, f"final={losses[-1]:.4f}")  # plot all three: one chart, one story