Performance, testing, and style ¶
Alfred Galichon (NYU) ¶
'math+econ+code' masterclass series: fundamentals of research in python ¶
With python code examples ¶
© 2018–2026 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.
If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass series. https://www.math-econ-code.org/
Learning objectives¶
Measure before optimizing, at three resolutions:
%timeitfor expressions,cProfilefor functions, andline_profilerfor individual lines.JIT-compile a genuinely sequential loop with Numba, and recognize when that is and is not the right tool.
Write a
pytestsuite that tests invariants and closed forms, not just remembered outputs; use parametrization and float-aware comparison.Enforce style with
ruffand catch type errors withmypybefore running the code.Take the
mec_numericalmodule offd03and make it something a colleague could depend on, which is what the remaining seven lectures will assume you can do.
References¶
[K] Knuth, D. E. (1974). "Structured Programming with go to Statements." Computing Surveys 6(4): the source of the remark in §2.
[Wi] Wilson, G. et al. (2017). "Good Enough Practices in Scientific Computing." PLOS Computational Biology 13(6).
[N] Numba documentation. https://numba.readthedocs.io/.
[P] pytest documentation. https://docs.pytest.org/.
[R] Ruff documentation. https://docs.astral.sh/ruff/.
[M] mypy documentation. https://mypy.readthedocs.io/.
1. From working code to reliable code¶
You can now write Python that works. The gap between code that works on your machine, today, in this notebook and code a colleague can run next year, on a different operating system, and get your numbers is software engineering, and it divides in two.
This lecture takes the first half: making code fast, correct, readable, and maintainable. Four practices, four tools, applied to numerical code of exactly the kind fd03 §8 packaged into a module.
| Pillar | Tools | Where the series raised it |
|---|---|---|
| Performance: fast enough for the problem | %timeit, cProfile, line_profiler, Numba |
the iteration counts of fd01 §12; the loops of fd02 §8 |
| Correctness: verifiable, regression-proof | pytest |
every "check passed" line since fd01 |
| Readability: conformant style | ruff |
fd02 §13 on PEP 8 |
| Maintainability: types that catch bugs early | mypy |
fd03 §11 on @dataclass annotations |
The second half is reproducibility: the practices that let someone else obtain your numbers at all. It needs the data, the acquired sources and the randomness that the next four lectures introduce, so fd07 assembles it once those debts exist. Two of its seven items are already in this lecture: test the code, and the environment discipline of fd03 §8.
A note on where this sits. It is deliberately early. Profiling and testing are habits, and a habit formed over seven remaining lectures is worth more than a checklist read at the end of the series.
2. Profile before you optimize¶
Premature optimization is the root of all evil.: Knuth (1974)
Intuition about where time goes is usually wrong. The loop worth optimizing is the one a profiler pointed at. Three tools, at increasing resolution.
%timeit runs an expression repeatedly and reports the mean and standard deviation per loop, auto-tuning the number of repetitions. %%timeit does the same for a whole cell.
import math
import time
import subprocess
import sys
import tempfile
from pathlib import Path
import numpy as np
OUT_DIR = Path("generated")
OUT_DIR.mkdir(exist_ok=True)
def portable(text):
"""Strip machine-specific path prefixes from tool output before printing.
Profilers and test runners report absolute paths. Those are useful to you
and meaningless -- occasionally revealing -- to a reader, so normalize them
on the way out. Temp first: on Windows it sits underneath the home
directory.
"""
for prefix, label in ((tempfile.gettempdir(), "<tmp>"), (str(Path.home()), "~")):
for form in (prefix, prefix.replace(chr(92), "/")):
text = text.replace(form, label)
return text
N = 100_000
def sum_squares_python(n):
return sum(k * k for k in range(n))
def sum_squares_numpy(n):
# int64 explicitly: NumPy's default integer is 32-bit on Windows, and k*k
# silently overflows past k = 46,340 -- fast, wrong, and no warning; fd08 section 3 returns to this.
return int((np.arange(n, dtype=np.int64) ** 2).sum())
assert sum_squares_python(N) == sum_squares_numpy(N)
print(f"both give {sum_squares_python(N):,}")
both give 333,328,333,350,000
t_py = %timeit -o -q sum_squares_python(N)
t_np = %timeit -o -q sum_squares_numpy(N)
print(f"pure Python : {t_py.average * 1e3:7.2f} ms +/- {t_py.stdev * 1e3:.2f}")
print(f"NumPy : {t_np.average * 1e3:7.2f} ms +/- {t_np.stdev * 1e3:.2f}")
print(f"speedup : {t_py.average / t_np.average:7.1f}x")
pure Python : 13.37 ms +/- 0.20 NumPy : 0.46 ms +/- 0.05 speedup : 29.3x
cProfile answers the next question: not how long, but where. Use it when a routine is slow and you do not know which call inside it is responsible.
Its report names files by absolute path, which is why every tool output in this lecture is printed through the portable helper defined above. Output you intend to publish should not carry your home directory: it is noise to a reader, it makes two runs on two machines look different when they are not, and it discloses more about your filesystem than you probably intend.
import cProfile
import pstats
from io import StringIO
def outer():
return sum_squares_python(50_000) + sum_squares_numpy(50_000)
profiler = cProfile.Profile()
profiler.enable()
for _ in range(20):
outer()
profiler.disable()
stream = StringIO()
pstats.Stats(profiler, stream=stream).sort_stats("cumulative").print_stats(6)
print(portable("\n".join(stream.getvalue().splitlines()[:14])))
1000248 function calls (1000247 primitive calls) in 0.300 seconds
Ordered by: cumulative time
List reduced from 48 to 6 due to restriction <6>
ncalls tottime percall cumtime percall filename:lineno(function)
3/2 0.000 0.000 0.286 0.143 ~\anaconda3\Lib\site-packages\IPython\core\interactiveshell.py:3541(run_code)
2 0.014 0.007 0.286 0.143 {built-in method builtins.exec}
20 0.000 0.000 0.283 0.014 <tmp>\ipykernel_9580\2432541108.py:28(sum_squares_python)
20 0.148 0.007 0.282 0.014 {built-in method builtins.sum}
20 0.000 0.000 0.271 0.014 <tmp>\ipykernel_9580\806927163.py:5(outer)
1000020 0.135 0.000 0.135 0.000 <tmp>\ipykernel_9580\2432541108.py:29(<genexpr>)
Each row gives the call count, the time spent inside that function excluding sub-calls (tottime), and the time including them (cumtime). Read cumtime from the top to find where time goes overall; read tottime to find which leaf function is doing the work.
3. Line-level profiling¶
cProfile stops at the function boundary. When a function is slow and you need to know which line, use line_profiler.
Two practical notes, and the second one matters.
- Profile a function defined in a module on disk, not in a notebook cell, so the profiler can find its source. We write one with
%%writefile, as infd03§8. - Use the
LineProfilerAPI, not the%lprunmagic. The magic sends its report to IPython's pager, which is a no-op undernbconvert, so a notebook using%lprunexecutes with no output at all and the section renders empty in the exported HTML. A tool that only works interactively is a tool that does not survive a clean re-run, which is this lecture's whole subject.
%%writefile generated/mec_perf.py
"""Routines profiled in fd04."""
import numpy as np
def solow_path(nbt, s=0.25, delta=0.05, alpha=0.30, k0=0.1):
"""Discrete-time Solow accumulation (fd03 section 9), written as a loop."""
k_t = np.empty(nbt)
k_t[0] = k0
for t in range(1, nbt):
y = k_t[t - 1] ** alpha
k_t[t] = s * y + (1 - delta) * k_t[t - 1]
return k_t
Overwriting generated/mec_perf.py
sys.path.insert(0, str(OUT_DIR)) if str(OUT_DIR) not in sys.path else None
from mec_perf import solow_path
from line_profiler import LineProfiler
profiler = LineProfiler(solow_path)
profiler.runcall(solow_path, 200_000)
report = StringIO()
profiler.print_stats(stream=report)
print(portable(report.getvalue()))
Timer unit: 1e-07 s
Total time: 0.361156 s
File: ~\Dropbox\AGResearch\courses\m-e-c\___research-tools\generated\mec_perf.py
Function: solow_path at line 5
Line # Hits Time Per Hit % Time Line Contents
==============================================================
5 def solow_path(nbt, s=0.25, delta=0.05, alpha=0.30, k0=0.1):
6 """Discrete-time Solow accumulation (fd03 section 9), written as a loop."""
7 1 311.0 311.0 0.0 k_t = np.empty(nbt)
8 1 61.0 61.0 0.0 k_t[0] = k0
9 200000 606665.0 3.0 16.8 for t in range(1, nbt):
10 199999 1117495.0 5.6 30.9 y = k_t[t - 1] ** alpha
11 199999 1886927.0 9.4 52.2 k_t[t] = s * y + (1 - delta) * k_t[t - 1]
12 1 100.0 100.0 0.0 return k_t
The % Time column is the answer. Roughly half the time is in the accumulation line and a third in the exponentiation, and the loop overhead itself, the for statement, accounts for a sixth. That tells you two things at once: there is no single dominant line to fix, and the cost is spread across interpreted operations rather than concentrated in one expensive call. When the profile looks like this, many cheap lines, each executed hundreds of thousands of times, the fix is not to rewrite a line but to remove the interpreter, which is §4.
4. Numba: compiling a loop the interpreter cannot vectorize¶
Most numerical work in Python is made fast by pushing loops down into compiled array code, which is fd08's subject. The Solow recursion is the case where that is not available: $k_{t+1}$ depends on $k_t$, so the loop is genuinely sequential and there is no array expression to write.
Numba JIT-compiles such a function to machine code on its first call, keyed on the argument types, and reuses the compiled version afterwards. The @njit decorator requests it.
from numba import njit
@njit
def solow_path_jit(nbt, s=0.25, delta=0.05, alpha=0.30, k0=0.1):
k_t = np.empty(nbt)
k_t[0] = k0
for t in range(1, nbt):
k_t[t] = s * k_t[t - 1] ** alpha + (1 - delta) * k_t[t - 1]
return k_t
t0 = time.perf_counter()
solow_path_jit(100) # first call: compile
print(f"first call (compile + tiny run): {(time.perf_counter() - t0) * 1e3:7.1f} ms")
T = 2_000_000
t0 = time.perf_counter(); k_py = solow_path(T); t_py = time.perf_counter() - t0
t0 = time.perf_counter(); k_jit = solow_path_jit(T); t_jit = time.perf_counter() - t0
print(f"\npure Python ({T:,} steps): {t_py * 1e3:8.1f} ms")
print(f"Numba JIT ({T:,} steps): {t_jit * 1e3:8.1f} ms")
print(f"speedup: {t_py / t_jit:8.1f}x")
gap = np.abs(k_py - k_jit).max()
print(f"\nmax |python - jit| = {gap:.2e} (tolerance 1e-12)")
assert gap < 1e-12
print("check passed: identical answers, and the check is not optional --")
print("a faster wrong answer is the worst possible outcome.")
first call (compile + tiny run): 610.6 ms
pure Python (2,000,000 steps): 1817.3 ms Numba JIT (2,000,000 steps): 140.2 ms speedup: 13.0x max |python - jit| = 0.00e+00 (tolerance 1e-12) check passed: identical answers, and the check is not optional -- a faster wrong answer is the worst possible outcome.
An order of magnitude or more is typical, and on tight numerical kernels Numba routinely matches hand-written C. The costs are real and worth stating: a compilation pause on the first call, a restricted subset of Python (NumPy and scalars, not pandas or arbitrary objects), and code that is harder to debug because the fast path is no longer Python.
Three rules. @njit is the idiomatic spelling of @jit(nopython=True); an unsupported operation raises TypingError rather than silently falling back. @njit(cache=True) persists the compiled code between sessions. And the equality check above is not ceremony: the entire risk of an optimization is that it changes the answer, so every optimization gets a test comparing it against the slow version it replaces. That is the bridge to §5.
5. Testing with pytest¶
A test verifies code with code. In research it does three things: it documents what a function should do, it catches regressions when you reorganize, and it lets you change code without fear.
The conventions: files named test_*.py, functions named test_*, plain assert.
What to test in research code is the part worth thinking about. Testing that a function returns the number it returned yesterday is nearly worthless: it locks in whatever it did, right or wrong. Test instead the things you can derive independently:
- closed forms the code should reproduce;
- invariants and symmetries: an output that must not depend on an irrelevant input;
- limit cases: zero, identity, the steady state;
- cross-checks against a second algorithm.
Every "check passed" in this series is one of those four. The suite below writes them down properly.
%%writefile generated/test_mec_perf.py
"""Tests for mec_perf.solow_path -- invariants, not remembered numbers."""
import numpy as np
import pytest
from mec_perf import solow_path
PARAMS = [(0.25, 0.05, 0.30), (0.20, 0.10, 0.30), (0.30, 0.05, 0.50)]
def steady_state(s, delta, alpha):
return (s / delta) ** (1 / (1 - alpha))
@pytest.mark.parametrize("s,delta,alpha", PARAMS)
def test_length(s, delta, alpha):
"""The path has exactly the requested number of periods."""
assert len(solow_path(50, s=s, delta=delta, alpha=alpha)) == 50
@pytest.mark.parametrize("s,delta,alpha", PARAMS)
def test_steady_state_is_fixed(s, delta, alpha):
"""Started AT the steady state, the path never leaves it."""
k_star = steady_state(s, delta, alpha)
path = solow_path(100, s=s, delta=delta, alpha=alpha, k0=k_star)
assert path == pytest.approx(k_star, rel=1e-12)
@pytest.mark.parametrize("s,delta,alpha", PARAMS)
def test_monotone_from_below(s, delta, alpha):
"""Started BELOW the steady state, the path increases and does not overshoot."""
k_star = steady_state(s, delta, alpha)
path = solow_path(200, s=s, delta=delta, alpha=alpha, k0=0.1 * k_star)
assert np.all(np.diff(path) > 0)
assert path[-1] < k_star
# 200 periods, NOT 400: after 400 the error reaches the rounding floor and the
# measured ratio is noise. See the cell below the test run.
RATE_PERIODS = 200
def measured_rate(s, delta, alpha, periods=RATE_PERIODS):
k_star = steady_state(s, delta, alpha)
path = solow_path(periods, s=s, delta=delta, alpha=alpha, k0=0.99 * k_star)
return (path[-1] - k_star) / (path[-2] - k_star)
@pytest.mark.parametrize("s,delta,alpha", PARAMS)
def test_convergence_rate(s, delta, alpha):
"""Near k*, the error contracts at 1 - delta(1 - alpha) [fd03 section 9]."""
assert measured_rate(s, delta, alpha) == pytest.approx(
1 - delta * (1 - alpha), rel=1e-5)
def test_saving_rate_moves_level_not_rate():
"""s relocates the steady state but not the convergence rate [fd03 Ex 4]."""
rates = [measured_rate(s, 0.05, 0.30) for s in (0.15, 0.25, 0.40)]
assert rates[0] == pytest.approx(rates[1], rel=1e-5)
assert rates[1] == pytest.approx(rates[2], rel=1e-5)
Overwriting generated/test_mec_perf.py
result = subprocess.run(
[sys.executable, "-m", "pytest", "test_mec_perf.py", "-v", "--tb=short",
"-p", "no:cacheprovider"],
capture_output=True, text=True, cwd=str(OUT_DIR))
print(portable(result.stdout[-2200:]))
assert result.returncode == 0, "the test suite failed"
print("suite passed.")
============================= test session starts ============================= platform win32 -- Python 3.12.3, pytest-8.4.2, pluggy-1.6.0 -- ~\anaconda3\python.exe rootdir: ~\Dropbox\AGResearch\courses\m-e-c\___research-tools\generated plugins: anyio-4.2.0 collecting ... collected 13 items test_mec_perf.py::test_length[0.25-0.05-0.3] PASSED [ 7%] test_mec_perf.py::test_length[0.2-0.1-0.3] PASSED [ 15%] test_mec_perf.py::test_length[0.3-0.05-0.5] PASSED [ 23%] test_mec_perf.py::test_steady_state_is_fixed[0.25-0.05-0.3] PASSED [ 30%] test_mec_perf.py::test_steady_state_is_fixed[0.2-0.1-0.3] PASSED [ 38%] test_mec_perf.py::test_steady_state_is_fixed[0.3-0.05-0.5] PASSED [ 46%] test_mec_perf.py::test_monotone_from_below[0.25-0.05-0.3] PASSED [ 53%] test_mec_perf.py::test_monotone_from_below[0.2-0.1-0.3] PASSED [ 61%] test_mec_perf.py::test_monotone_from_below[0.3-0.05-0.5] PASSED [ 69%] test_mec_perf.py::test_convergence_rate[0.25-0.05-0.3] PASSED [ 76%] test_mec_perf.py::test_convergence_rate[0.2-0.1-0.3] PASSED [ 84%] test_mec_perf.py::test_convergence_rate[0.3-0.05-0.5] PASSED [ 92%] test_mec_perf.py::test_saving_rate_moves_level_not_rate PASSED [100%] ============================= 13 passed in 0.15s ============================== suite passed.
The first version of that suite failed, and the reason is worth the detour. I originally measured the convergence rate after 400 periods, reasoning that more iterations means closer to the asymptotic regime. Two tests failed. Here is why.
k_star_demo = (0.20 / 0.10) ** (1 / (1 - 0.30)) # s=0.20, delta=0.10, alpha=0.30
predicted = 1 - 0.10 * (1 - 0.30)
print(f"predicted contraction rate: {predicted}\n")
print(f"{'periods':>9}{'error':>13}{'relative':>12}{'measured ratio':>18}{'gap':>10}")
print("-" * 62)
for periods in (30, 60, 100, 200, 400):
path = solow_path(periods, s=0.20, delta=0.10, alpha=0.30, k0=0.99 * k_star_demo)
e_prev, e_last = path[-2] - k_star_demo, path[-1] - k_star_demo
ratio = e_last / e_prev
print(f"{periods:>9}{abs(e_prev):>13.2e}{abs(e_prev) / k_star_demo:>12.1e}"
f"{ratio:>18.9f}{abs(ratio - predicted):>10.1e}")
print("-" * 62)
print(f"at 400 periods the error is ~1e-14 -- at the rounding floor -- and the")
print(f"ratio collapses to {16/17:.9f} = 16/17, a quotient of two numbers that are")
print(f"pure floating-point noise. The rate is unmeasurable there.")
predicted contraction rate: 0.93
periods error relative measured ratio gap
--------------------------------------------------------------
30 3.53e-03 1.3e-03 0.930013793 1.4e-05
60 4.01e-04 1.5e-04 0.930001563 1.6e-06
100 2.20e-05 8.2e-06 0.930000086 8.6e-08
200 1.55e-08 5.8e-09 0.929999992 8.3e-09
400 7.55e-15 2.8e-15 0.941176471 1.1e-02
--------------------------------------------------------------
at 400 periods the error is ~1e-14 -- at the rounding floor -- and the
ratio collapses to 0.941176471 = 16/17, a quotient of two numbers that are
pure floating-point noise. The rate is unmeasurable there.
This is the trap fd01 §12 identified, and it caught the
author of both while writing this lecture's test suite. A rate can only be measured while the
quantity whose decay defines it is still above the noise floor; iterating further eventually
destroys the measurement rather than improving it.
Two lessons follow. First, a test asserting an asymptotic property must specify the window in which
the assertion is meaningful, which is why RATE_PERIODS = 200 is a named constant with a comment
rather than a number buried in three tests. Second, and more useful: the failing test was correct and the reasoning was wrong. The suite caught an error in the test rather than in the code, which is
what a test suite is for, and a good reason to be suspicious of the reflex to loosen a tolerance
until it passes.
Note how the suite is invoked: subprocess.run([sys.executable, "-m", "pytest", ...]) rather than the !pytest shell escape. sys.executable is this kernel's interpreter, so the tests necessarily run against the same environment as the notebook. A bare !pytest runs whatever the shell's PATH finds first, which on a machine with several environments is frequently not the one you are working in: the interpreter-versus-kernel distinction of fd01 §4, in a particularly consequential form.
Look at what the tests assert. Not one of them hard-codes a number the code once produced. test_steady_state_is_fixed asserts a fixed point; test_monotone_from_below asserts monotonicity and no overshoot; test_convergence_rate asserts the rate derived in fd03 §9; and the last asserts the $s$-independence proved in fd03, Exercise 4. Each would catch a wrong implementation that nevertheless ran. That is the difference between a test suite and a snapshot.
Three features useful: @pytest.mark.parametrize turns one test into many, each reported separately; @pytest.fixture injects shared setup; and pytest.approx compares floats with a tolerance, because fd01 §5.
6. Style: ruff¶
Consistent style is not aesthetics. It removes friction when reading code: a collaborator's, or your own from a year ago. ruff format formats in place and ruff check lints, both written in Rust and fast enough to run on every save.
%%writefile generated/ugly.py
import numpy as np
import os
def f( x,y ):
z = x +y
return z* 2
class Thing :
def __init__( self,name):
self.name=name
def greet(self ):
print( "hi "+self.name )
unused = "never referenced again"
Overwriting generated/ugly.py
lint = subprocess.run([sys.executable, "-m", "ruff", "check", str(OUT_DIR / "ugly.py")],
capture_output=True, text=True)
print("ruff check:\n" + (lint.stdout or lint.stderr)[:800])
subprocess.run([sys.executable, "-m", "ruff", "format", str(OUT_DIR / "ugly.py")],
capture_output=True, text=True)
print("\nafter ruff format:\n")
print((OUT_DIR / "ugly.py").read_text())
ruff check: I001 [*] Import block is un-sorted or un-formatted --> generated\ugly.py:1:1 | 1 | / import numpy as np 2 | | import os | |_________^ 3 | def f( x,y ): 4 | z = x +y | help: Organize imports | - import numpy as np 1 | import os 2 + 3 + import numpy as np 4 + 5 + 6 | def
after ruff format:
import numpy as np
import os
def f(x, y):
z = x + y
return z * 2
class Thing:
def __init__(self, name):
self.name = name
def greet(self):
print("hi " + self.name)
unused = "never referenced again"
ruff check found the unused import that ruff format will not touch: the two do different jobs. Formatting is mechanical and safe; linting reports things that might be bugs and require a decision. Run the formatter automatically; read the linter.
7. Type hints and mypy¶
Annotations document the expected types. Python ignores them at runtime, but mypy reads them and reports mismatches before the code runs.
%%writefile generated/newton_typed.py
"""A typed Newton implementation, with one deliberate error."""
from typing import Callable
def newton(
f: Callable[[float], float],
df: Callable[[float], float],
x0: float,
tol: float = 1e-10,
maxiter: int = 100,
) -> tuple[float, int]:
"""Solve f(x) = 0 by Newton's method. Returns (root, iterations)."""
x = x0
for k in range(maxiter):
fx = f(x)
if abs(fx) < tol:
return x, k
x = x - fx / df(x)
return x, maxiter
# a string where a float is expected -- mypy should catch this without running it
result: tuple[float, int] = newton(lambda x: x**2 - 2, lambda x: 2 * x,
x0="this should be a float")
Overwriting generated/newton_typed.py
check = subprocess.run([sys.executable, "-m", "mypy", str(OUT_DIR / "newton_typed.py")],
capture_output=True, text=True)
print(check.stdout or check.stderr)
assert "error:" in check.stdout, "mypy was expected to find the deliberate error"
print("check passed: the type error was found without executing the file.")
generated\newton_typed.py:24: error: Argument "x0" to "newton" has incompatible type "str"; expected "float" [arg-type] Found 1 error in 1 file (checked 1 source file) check passed: the type error was found without executing the file.
A TypeError that would have surfaced at runtime, possibly only on a rarely-exercised code path, became a located diagnostic at edit time.
A practical rule: annotate function signatures and class attributes; leave local variables to inference. X | None for optional, Callable[[float], float] for a function argument, Any as the escape hatch. And note that hints help readers even with no checker in sight: def estimate(X: np.ndarray, y: np.ndarray, *, ridge: float = 0.0) -> np.ndarray tells you compactly what the function consumes and returns.
8. Summary¶
Measure, then optimize.
%timeitfor an expression,cProfilefor a call tree,line_profilerfor a line, and use its API rather than%lprun, which produces no output outside an interactive session.Numba is for the loop you cannot vectorize. It bought better than an order of magnitude on the Solow recursion, which is genuinely sequential. Every optimization is checked against the implementation it replaces, because a faster wrong answer is an undesirable outcome.
Test invariants, not remembered outputs. The suite in §5 asserts a fixed point, monotonicity, a convergence rate derived in
fd03§9, and an $s$-independence proved infd03Exercise 4. A test that hard-codes yesterday's number locks in whatever the code did, correct or not.A failing test is evidence, not an obstacle. The one in §5 was right and the reasoning behind it was wrong: the rate had been measured below the rounding floor. Loosening a tolerance until a test passes discards exactly the information you asked for.
ruffandmypycost nothing per run and pay for themselves at edit time. The formatter is mechanical and safe; the linter reports things that require a decision; the type checker finds errors on code paths you have not executed.What is deferred. The other half of shipping code is reproducibility: environments, seeds, data manifests, provenance.
fd07§8 assembles it, once there is data and randomness for it to apply to.
9. Exercises¶
Worked solutions are in §11.
Exercise 1: Find the hot line, then remove it. fd01 §10 priced a cash-flow stream. Write a module generated/mec_pv.py containing
as an explicit Python loop that carries the discount factor forward by dividing, rather than recomputing a power each step.
Profile it with LineProfiler at $T = 200{,}000$ and identify the line taking the most time. Then write a vectorized version and an @njit version, and check all three against the closed form for a level perpetuity, $PV = c(1+r)/r$.
Which is faster, vectorizing or Numba, and why? Then one more measurement: run the same loop over a Python list instead of a NumPy array, and explain the result.
# your answer here
Exercise 2: When Numba does not help. @njit gave a large speedup on the Solow loop. Find two cases where it does not.
- Apply
@njitto a function that is already a single NumPy expression: saylambda x: (x ** 2).sum()on an array of $10^7$, and compare. - Apply it to a function doing string or
dictwork, and report what happens.
Explain both results in terms of what Numba actually removes. Then state the rule you would give a colleague for deciding whether to use it.
# your answer here
Exercise 3: A test that would have caught a real bug. The Solow steady state solves $s k^{\alpha} - \delta k = 0$ for $k > 0$. With $s = 0.25$, $\delta = 0.05$, $\alpha = 0.30$ the answer is $k^{\star} = (s/\delta)^{1/(1-\alpha)} \approx 9.9662$.
Call newton from fd03's mec_numerical on this equation starting from $k_0 = 1$. It reports converged=True. Inspect what it actually returned, and explain in one sentence what happened at the first step.
Now write a pytest file with two tests: one comparing the root to the closed form with pytest.approx, one asserting that the returned root is a real number, and run it against three solvers: newton from $k_0=1$; a damped Newton that halves the step until the iterate stays positive; and bisect on the bracket $[10^{-6}, 100]$.
Report which tests each solver passes. The point to extract: neither test alone is sufficient, and each catches what the other misses.
# your answer here
Exercise 4: The cost of the wrong container. fd02 §6 said membership testing is $O(n)$ in a list and $O(1)$ in a set. Measure it.
For $n \in \{10^2, 10^3, 10^4, 10^5, 10^6\}$, build list(range(n)) and set(range(n)) and time a membership test for an element that is absent: the worst case, and the one that matters, since a lookup that fails must scan the whole list.
Tabulate nanoseconds per lookup for both, and the ratio. Confirm the two growth rates from the numbers rather than from the claim, and say what the measurement implies for a loop that filters one collection against another.
# your answer here
Exercise 5: How badly does brute force scale? fd02 §12 solved a $3\times4$ assignment problem by enumerating all $24$ matchings, and remarked that this is "astronomical at realistic sizes." Put a number on it.
For $n = 4,\dots,10$, build an $n\times n$ surplus matrix from any deterministic rule, find the optimal assignment by enumerating all $n!$ permutations, and time it. Tabulate $n$, $n!$, and the time per permutation.
Then extrapolate to $n = 15$, $20$ and $25$ using the measured cost per permutation, and report the wall-clock in units a reader can feel. Finally: what does this tell you about the relationship between the statement of an optimization problem and its solvability, and why is fd09 §5 going to matter?
# your answer here
10. Further directions¶
You can now write Python that is fast enough, tested, styled, and type-checked. What is still missing is everything that makes it someone else's, and that begins with version control.
fd05 puts the code in a repository: the commit graph, branching and merging, the pull-request workflow, and the research compendium that packages code, environment and data as one object. The tests/ folder in that compendium is the suite you just wrote.
From there, fd06 and fd07 bring in data, tables, and the web sources they come from, and fd07 §8 closes the second half of this lecture by assembling reproducibility as seven concrete practices. Two of them are already yours.
11. Solutions to the exercises¶
Solution to Exercise 1: Find the hot line, then remove it¶
%%writefile generated/mec_pv.py
"""Present value of a cash-flow stream, as a loop (fd01 section 10)."""
def pv_loop(c_t, r):
"""PV = sum_t c_t (1+r)^{-t}, carrying the discount factor forward."""
total = 0.0
disc = 1.0
for t in range(len(c_t)):
total += c_t[t] * disc
disc /= 1.0 + r
return total
Overwriting generated/mec_pv.py
import importlib
import mec_pv
importlib.reload(mec_pv)
from mec_pv import pv_loop
T_PV, R, C = 200_000, 0.03, 100.0
c_t = np.full(T_PV, C)
exact = C * (1.0 + R) / R # level perpetuity, fd01 section 11
prof = LineProfiler(pv_loop)
prof.runcall(pv_loop, c_t, R)
rep = StringIO()
prof.print_stats(stream=rep)
print(portable(rep.getvalue()))
Timer unit: 1e-07 s
Total time: 0.268698 s
File: ~\Dropbox\AGResearch\courses\m-e-c\___research-tools\generated\mec_pv.py
Function: pv_loop at line 4
Line # Hits Time Per Hit % Time Line Contents
==============================================================
4 def pv_loop(c_t, r):
5 """PV = sum_t c_t (1+r)^{-t}, carrying the discount factor forward."""
6 1 15.0 15.0 0.0 total = 0.0
7 1 3.0 3.0 0.0 disc = 1.0
8 200001 591476.0 3.0 22.0 for t in range(len(c_t)):
9 200000 1352845.0 6.8 50.3 total += c_t[t] * disc
10 200000 742535.0 3.7 27.6 disc /= 1.0 + r
11 1 108.0 108.0 0.0 return total
Half the time is in the accumulation and a quarter in the division, with the for statement itself taking the rest. Now remove the interpreter, two ways.
def pv_vec(c_t, r):
return float((c_t * (1.0 + r) ** -np.arange(len(c_t))).sum())
@njit
def pv_jit(c_t, r):
total = 0.0
disc = 1.0
for t in range(len(c_t)):
total += c_t[t] * disc
disc /= 1.0 + r
return total
pv_jit(c_t[:10], R) # compile
print(f"{'method':<14}{'time (ms)':>12}{'PV':>18}{'|PV - c(1+r)/r|':>20}{'speedup':>10}")
print("-" * 74)
timings = {}
for name, fn in (("loop", pv_loop), ("vectorized", pv_vec), ("numba", pv_jit)):
t0 = time.perf_counter()
v = fn(c_t, R)
dt = time.perf_counter() - t0
timings[name] = dt
print(f"{name:<14}{dt * 1e3:>12.2f}{v:>18.10f}{abs(v - exact):>20.2e}"
f"{timings['loop'] / dt:>10.1f}x")
print("-" * 74)
print(f"closed form c(1+r)/r = {exact:.10f}")
method time (ms) PV |PV - c(1+r)/r| speedup --------------------------------------------------------------------------
loop 99.69 3433.3333333333 1.36e-11 1.0x vectorized 6.64 3433.3333333333 3.64e-12 15.0x numba 0.65 3433.3333333333 1.36e-11 153.6x -------------------------------------------------------------------------- closed form c(1+r)/r = 3433.3333333333
Numba wins here, and by a wide margin over vectorizing, which is the opposite of the ordering §4 would lead you to expect, and the reason is worth having.
pv_vec must materialize what the loop only ever holds one number at a time: an array of $T$ discount factors, then an array of $T$ products, then a reduction. That is three passes over 1.6 MB of memory for a computation whose entire state is two scalars. The @njit loop allocates nothing, touches c_t once, and keeps total and disc in registers. When an algorithm is genuinely sequential and its state is $O(1)$, the compiled loop is not merely competitive with the array form: it is the better algorithm, and vectorizing it is a pessimization dressed as an optimization.
And the accuracy differs, in the direction you would not guess. The two loops carry the discount factor forward by repeated division, so $T$ roundings accumulate; pv_vec computes each $(1+r)^{-t}$ independently in one operation. The vectorized error is roughly a quarter of the loops'. The fast version is the less accurate one. Both are far inside any tolerance that matters here, but the trade is real and it is the sort of thing to know you are making.
lst = [C] * T_PV # identical values, Python list
best = {}
for name, obj in (("ndarray", c_t), ("list", lst)):
reps = []
for _ in range(3):
t0 = time.perf_counter(); pv_loop(obj, R); reps.append(time.perf_counter() - t0)
best[name] = min(reps)
print(f"pv_loop over a {name:<9}{best[name] * 1e3:8.2f} ms")
print(f"\nthe SAME loop is {best['ndarray'] / best['list']:.1f}x slower over the array")
pv_loop over a ndarray 95.66 ms pv_loop over a list 34.47 ms the SAME loop is 2.8x slower over the array
The last measurement is the trap. The identical Python loop runs roughly three times faster over a plain list than over the NumPy array, so a colleague who "speeds up" a loop by first putting the data in an array has made it slower.
The reason is boxing. lst[t] returns a reference to a Python float object that already exists; c_t[t] must construct a new Python float wrapping the raw 8 bytes at that offset, allocate it, and free it next iteration. NumPy's speed comes entirely from operations that never leave compiled code; index it one element at a time from Python and you pay for the packing without ever using it. The rule that follows is sharp: an array is either processed whole or not put in an array at all.
Solution to Exercise 2: When Numba does not help¶
big = np.linspace(-3.0, 3.0, 10_000_000)
def sum_sq_numpy(x):
return (x ** 2).sum()
@njit
def sum_sq_jit(x):
return (x ** 2).sum()
sum_sq_jit(big[:10]) # compile
t0 = time.perf_counter(); a = sum_sq_numpy(big); t_np = time.perf_counter() - t0
t0 = time.perf_counter(); b = sum_sq_jit(big); t_jit = time.perf_counter() - t0
print(f"{'already-vectorized expression':<32}")
print(f" NumPy : {t_np * 1e3:8.1f} ms")
print(f" Numba : {t_jit * 1e3:8.1f} ms ({t_np / t_jit:.2f}x)")
print(f" agree : {abs(a - b) / abs(a):.1e} relative")
already-vectorized expression NumPy : 44.1 ms Numba : 53.5 ms (0.82x) agree : 3.7e-16 relative
import re
import warnings
# 2. an untyped dict built from an empty literal
@njit
def count_words_jit(words):
counts = {}
for w in words:
counts[w] = counts.get(w, 0) + 1
return counts
def plain(text): # Numba colours its messages; strip the codes
return re.sub(r"\x1b\[[0-9;]*m", "", text)
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
try:
count_words_jit(["alpha", "beta", "alpha"])
except Exception as err:
print(f"{type(err).__name__}: {plain(str(err)).splitlines()[0][:90]}")
for w in caught:
print(f" warning: {w.category.__name__}: "
f"{plain(str(w.message)).splitlines()[0][:72]}")
TypingError: Failed in nopython mode pipeline (step: nopython frontend) warning: NumbaTypeSafetyWarning: unsafe cast from unicode_type to undefined. Precision may be lost. warning: NumbaTypeSafetyWarning: unsafe cast from int64 to undefined. Precision may be lost.
But do not read that as "Numba cannot do strings or dictionaries." It can do both, and the distinction matters when you are deciding whether to use it.
from numba.typed import Dict
from numba.core import types
@njit
def shout(s): # unicode strings: supported
return s.upper() + "!"
@njit
def count_typed(words, counts): # a dict with declared types
for w in words:
counts[w] = counts.get(w, 0) + 1
return counts
typed = Dict.empty(key_type=types.unicode_type, value_type=types.int64)
print(f"unicode in nopython mode : {shout('hello')}")
print(f"typed dict in nopython : {dict(count_typed(['alpha', 'beta', 'alpha'], typed))}")
print("\nso the failure above is about type *inference*, not about strings or dicts.")
unicode in nopython mode : HELLO!
typed dict in nopython : {'alpha': 2, 'beta': 1}
so the failure above is about type *inference*, not about strings or dicts.
@njit
def build_inferable(): # an empty {} literal -- but typed
d = {} # by the assignment on the next line
d["alpha"] = 1
d["beta"] = 2
return d
print(f"empty literal, types inferable : {dict(build_inferable())}")
print("even `{}` compiles when later use determines the key and value types;")
print("what defeats it above is `.get` on a dict whose types are still unknown.")
empty literal, types inferable : {'alpha': 1, 'beta': 2}
even `{}` compiles when later use determines the key and value types;
what defeats it above is `.get` on a dict whose types are still unknown.
Why neither case improves performance. Numba removes the Python interpreter from a loop. In the first case there is no Python loop to remove: (x ** 2).sum() already executes entirely inside NumPy's compiled code, so there is nothing left for Numba to take away. It re-implements the same compiled work and does it worse, about a third slower here, because NumPy's reductions are hand-tuned and pairwise-summed, which a naive compiled loop is not. In the second, the obstacle is type inference, not the data types themselves. counts = {} gives Numba an empty dictionary with nothing to infer key and value types from, and counts.get(w, 0) on that object leaves them undetermined, so compilation fails with TypingError. The NumbaTypeSafetyWarning lines printed underneath are part of that failure rather than a separate problem: the inferencer reporting casts it could not resolve on its way to giving up. As the second cell shows, Numba compiles Unicode strings happily, and compiles the same counting loop once the dictionary arrives as a numba.typed.Dict with declared key and value types. What it refuses is the dynamically typed Python container, and it refuses rather than quietly reverting to interpreted execution: a feature, since the old "object mode" fallback produced code that compiled, ran, and was no faster, an undesirable compromise.
Practical rule. Reach for Numba when a profiler shows time inside an explicit Python loop over scalars that you cannot express as array operations, typically because each iteration depends on the last, as in the Solow recursion, or because the array form would allocate memory the loop does not need, as in Exercise 1. If the hot code is already a NumPy expression, Numba has nothing to offer; if it relies on containers whose types it cannot pin down, an empty dictionary whose key or value types are never resolved, a list of mixed types, a DataFrame, an arbitrary Python object, it may not compile, and the fix is usually to make the types explicit rather than to abandon Numba. Note may: as the cell below shows, {} is fine whenever later use determines its types. See the Numba reference on supported Python features [N].
Solution to Exercise 3: A test that would have caught a real bug¶
%%writefile generated/mec_solvers.py
"""Three ways to solve s k^alpha - delta k = 0 for k > 0.
`newton` and `bisect` are reproduced from fd03's mec_numerical so this file
stands alone; `newton_damped` is the guarded variant of Exercise 3.
"""
def newton(f, df, x0, tol=1e-10, maxiter=100):
x = x0
for k in range(maxiter):
fx = f(x)
if abs(fx) < tol:
return x, k, True
x = x - fx / df(x)
return x, maxiter, False
def newton_damped(f, df, x0, in_domain, tol=1e-10, maxiter=100):
"""Newton, halving the step until the next iterate stays in the domain."""
x = float(x0)
for k in range(maxiter):
fx = f(x)
if abs(fx) < tol:
return x, k, True
step = fx / df(x)
for _ in range(60):
if in_domain(x - step):
break
step *= 0.5
else:
return x, k, False
x = x - step
return x, maxiter, False
def bisect(f, a, b, tol=1e-10, maxiter=200):
fa = f(a)
if fa * f(b) > 0:
raise ValueError("f(a) and f(b) must have opposite signs")
for k in range(maxiter):
m = 0.5 * (a + b)
fm = f(m)
if abs(fm) < tol or 0.5 * (b - a) < tol:
return m, k, True
if fa * fm < 0:
b = m
else:
a, fa = m, fm
return m, maxiter, False
Overwriting generated/mec_solvers.py
from mec_solvers import newton # fd03's newton, verbatim
S, DELTA, ALPHA = 0.25, 0.05, 0.30
K_STAR = (S / DELTA) ** (1 / (1 - ALPHA))
f = lambda k: S * k ** ALPHA - DELTA * k
df = lambda k: S * ALPHA * k ** (ALPHA - 1) - DELTA
root, iters, converged = newton(f, df, 1.0)
print(f"k* (closed form) = {K_STAR!r}")
print(f"newton from k0=1 -> {root!r}")
print(f" converged={converged}, {iters} iterations, "
f"type={type(root).__name__}")
k* (closed form) = 9.966176578193442
newton from k0=1 -> (9.966176578162179-6.431023343521879e-12j)
converged=True, 6 iterations, type=complex
What happened at the first step. $f(1) = 0.20$ and $f'(1) = 0.025$, so Newton's step is $1 - 0.20/0.025 = -7$: the iterate leaves the domain in one move. Python then evaluates $(-7)^{0.3}$, and rather than raising, it returns a complex number: the principal branch of $z^{0.3}$. From there the iteration proceeds in complex arithmetic, abs() of a complex number is its modulus so the convergence test still works, and the solver reports success while returning a root with an imaginary part of about $-6\times10^{-12}$.
Nothing raises. The failure surfaces later and elsewhere: at the first float(root), or the first comparison root > 0, both of which raise TypeError a long way from the cause.
%%writefile generated/test_solvers.py
"""Two invariants, three solvers -- and neither test is sufficient alone."""
import pytest
from mec_solvers import bisect, newton, newton_damped
S, DELTA, ALPHA = 0.25, 0.05, 0.30
K_STAR = (S / DELTA) ** (1 / (1 - ALPHA))
def f(k):
return S * k ** ALPHA - DELTA * k
def df(k):
return S * ALPHA * k ** (ALPHA - 1) - DELTA
SOLVERS = {
"newton": lambda: newton(f, df, 1.0),
"damped": lambda: newton_damped(f, df, 1.0, lambda k: k > 0),
"bisect": lambda: bisect(f, 1e-6, 100.0),
}
@pytest.mark.parametrize("name", list(SOLVERS))
def test_matches_the_closed_form(name):
root, _, _ = SOLVERS[name]()
assert root == pytest.approx(K_STAR, rel=1e-6)
@pytest.mark.parametrize("name", list(SOLVERS))
def test_root_is_a_real_number(name):
root, _, _ = SOLVERS[name]()
assert isinstance(root, float)
Overwriting generated/test_solvers.py
result = subprocess.run(
[sys.executable, "-m", "pytest", "test_solvers.py", "-v", "--tb=no",
"-p", "no:cacheprovider"],
capture_output=True, text=True, cwd=str(OUT_DIR))
for line in result.stdout.splitlines():
if "::" in line or line.startswith("=") or line.startswith("FAILED"):
print(line)
test_solvers.py::test_matches_the_closed_form[newton] PASSED [ 16%] test_solvers.py::test_matches_the_closed_form[damped] FAILED [ 33%] test_solvers.py::test_matches_the_closed_form[bisect] PASSED [ 50%] test_solvers.py::test_root_is_a_real_number[newton] FAILED [ 66%] test_solvers.py::test_root_is_a_real_number[damped] PASSED [ 83%] test_solvers.py::test_root_is_a_real_number[bisect] PASSED [100%] FAILED test_solvers.py::test_matches_the_closed_form[damped] - assert 3.777778311260066e-32 == 9.966176578193442 ± 1.0e-05 FAILED test_solvers.py::test_root_is_a_real_number[newton] - assert False
from mec_solvers import bisect, newton_damped
print(f"{'solver':<24}{'root':>36}{'real?':>8}{'= k*?':>8}")
print("-" * 76)
for name, call in (("newton, k0=1", lambda: newton(f, df, 1.0)),
("newton damped, k0=1", lambda: newton_damped(f, df, 1.0, lambda k: k > 0)),
("bisect on [1e-6, 100]", lambda: bisect(f, 1e-6, 100.0))):
r, _, _ = call()
real = isinstance(r, float)
near = str(abs(r - K_STAR) < 1e-6) if real else "n/a"
print(f"{name:<24}{r!r:>36}{real!s:>8}{near:>8}")
print("-" * 76)
print(f"k* = {K_STAR!r}")
solver root real? = k*? ---------------------------------------------------------------------------- newton, k0=1 (9.966176578162179-6.431023343521879e-12j) False n/a newton damped, k0=1 3.777778311260066e-32 True False bisect on [1e-6, 100] 9.966176577593039 True True ---------------------------------------------------------------------------- k* = 9.966176578193442
Read the table by columns, because that is the lesson.
| solver | returns | test_matches_the_closed_form |
test_root_is_a_real_number |
|---|---|---|---|
newton from $k_0=1$ |
a complex number near $k^{\star}$ | passes | fails |
| damped Newton | $\approx 4\times10^{-32}$ | fails | passes |
bisect on a bracket |
$9.9661765776$ | passes | passes |
The value test passes on the complex root, because pytest.approx compares complex numbers by modulus and the imaginary part is $10^{-12}$. A suite consisting only of "does it get the right number" certifies a function that is silently returning the wrong type, and the eventual TypeError gets blamed on whatever code first tried to use the result.
The type test, on its own, is equally insufficient. The damped Newton never leaves the positive reals, and walks to $k = 0$, which is a perfectly good root of $sk^{\alpha} - \delta k$ and economically meaningless: the degenerate no-capital steady state. It returns a real, positive, correct-by-its-own-lights answer to the wrong question.
Only bisect passes both, and not by luck. It cannot fail either test by construction: it returns the midpoint of a bracket it never leaves, so the answer is real and inside $[10^{-6}, 100]$ by definition, and the sign change guarantees a root is in there. The lesson is not "write more tests" but the older one, when a method can be chosen that makes a failure mode impossible rather than detected, choose it, and let the tests guard the invariants that remain.
Solution to Exercise 4: The cost of the wrong container¶
import timeit
print(f"{'n':>10}{'list (ns)':>14}{'set (ns)':>12}{'ratio':>12}{'ns per element':>17}")
print("-" * 65)
for n in (100, 1_000, 10_000, 100_000, 1_000_000):
lst, st = list(range(n)), set(range(n))
absent = -1
reps = max(3, min(2000, 20_000_000 // n))
t_l = min(timeit.repeat(lambda: absent in lst, number=reps, repeat=3)) / reps
t_s = min(timeit.repeat(lambda: absent in st, number=reps, repeat=3)) / reps
print(f"{n:>10,}{t_l * 1e9:>14,.0f}{t_s * 1e9:>12,.0f}{t_l / t_s:>12,.0f}"
f"{t_l / n * 1e9:>17.2f}")
n list (ns) set (ns) ratio ns per element
-----------------------------------------------------------------
100 1,272 100 13 12.72
1,000 12,357 100 124 12.36
10,000 124,977 100 1,250 12.50
100,000 1,298,923 101 12,861 12.99
1,000,000 12,555,990 110 114,145 12.56
Both growth rates are visible in the numbers rather than assumed. The list column multiplies by ten each time $n$ does, that is $O(n)$, and the last column, time divided by $n$, is flat at about ten nanoseconds per element scanned, which is the constant hiding inside the $O$. The set column does not move at all: seventy to a hundred nanoseconds whether the set holds a hundred elements or a million, because hashing the key and probing one bucket does not depend on how many other buckets exist.
At $n = 10^6$ the ratio is about $10^5$. That is not a micro-optimization; it is the difference between a script that returns and one that does not.
What it implies for filtering one collection against another. The idiom [x for x in a if x not in b] costs $O(|a|\,|b|)$ when b is a list and $O(|a|)$ when b is a set, and the fix is one call, b = set(b), paid once. This is the most common accidental quadratic in research code, and it is invisible during development because it only bites at full data size. fd02 §6 gave the rule; this is why it is worth obeying without thinking about it.
Solution to Exercise 5: How badly does brute force scale?¶
import math
from itertools import permutations
def surplus(i, j):
return ((i + 1) * (j + 2)) % 7 + 1.0 # deterministic, no seed needed
print(f"{'n':>4}{'n!':>16}{'time (s)':>12}{'ns / permutation':>20}")
print("-" * 52)
cost = None
for n in range(4, 11):
Phi = [[surplus(i, j) for j in range(n)] for i in range(n)]
t0 = time.perf_counter()
best = max(sum(Phi[i][p[i]] for i in range(n)) for p in permutations(range(n)))
dt = time.perf_counter() - t0
cost = dt / math.factorial(n)
print(f"{n:>4}{math.factorial(n):>16,}{dt:>12.4f}{cost * 1e9:>20.1f}")
print("-" * 52)
def readable(seconds):
for unit, size in (("seconds", 1), ("hours", 3600), ("days", 86_400),
("years", 3.15576e7)):
if seconds / size < 400 or unit == "years":
return f"{seconds / size:,.1f} {unit}"
print(f"\nextrapolating at {cost * 1e9:.1f} ns per permutation:")
for n in (12, 15, 20, 25):
print(f" n = {n:<3}{math.factorial(n):>28,} permutations -> "
f"{readable(cost * math.factorial(n)):>22}")
n n! time (s) ns / permutation ---------------------------------------------------- 4 24 0.0000 1970.8 5 120 0.0002 1274.2 6 720 0.0010 1357.2 7 5,040 0.0075 1491.3 8 40,320 0.0650 1613.1
9 362,880 0.6354 1751.1
10 3,628,800 6.6263 1826.0 ---------------------------------------------------- extrapolating at 1826.0 ns per permutation: n = 12 479,001,600 permutations -> 0.2 hours n = 15 1,307,674,368,000 permutations -> 27.6 days n = 20 2,432,902,008,176,640,000 permutations -> 140,775.0 years n = 25 15,511,210,043,330,985,984,000,000 permutations -> 897,524,906,577.3 years
The cost per permutation barely moves. Ignore the $n=4$ row, which is timing noise on a loop of twenty-four; from there it drifts up by about a third across six orders of magnitude, because each evaluation sums $n$ terms. Essentially all of the growth is $n!$ itself, and $n!$ is not a function that rewards patience. A $20\times20$ assignment problem is small by any economic standard: twenty workers, twenty firms, a table that fits on one page. Enumerating it takes on the order of a hundred thousand years. At $n = 25$ it is some fifty times the age of the universe.
The statement of a problem and its solvability are different questions. Nothing about "choose the matching that maximizes total surplus" hints that the obvious algorithm is unusable at the size of a seminar room. The combinatorial explosion lives in the feasible set, not in the objective, and it is invisible in the way the problem is written down.
What rescues it is a change of representation. fd09 §5 writes the same problem as a linear program over the matching masses $\mu_{xy}$ rather than over permutations: the feasible set becomes a polytope cut out by $2n$ linear constraints, its vertices turn out to be exactly the permutations, and an interior-point method finds the optimum in time polynomial in $n$. A $20\times20$ instance is then milliseconds.
The change of representation pays a second dividend, which is the one this series cares about. That linear program has a dual, and fd09 §6 reads its multipliers as wages: a price system that implements the optimal assignment without anyone enumerating anything. The reason enumeration is hopeless and the reason prices exist turn out to be the same fact about the problem.