ulam

A general-purpose Monte Carlo simulation engine that scales from a thousand trials to billions - on a laptop.

Named for Stanisław Ulam, who co-invented the Monte Carlo method at Los Alamos in 1946 while wondering about the odds of winning a game of solitaire. The idea he and John von Neumann built out of that question - estimate what you can't solve exactly by sampling it, over and over, and let the law of large numbers do the rest - still works exactly the same way today, for anything from a coin flip to a billion-scenario risk model.

ulam is a from-scratch implementation of that idea, built to impose exactly one requirement on whatever you're simulating: give it a function shaped model(rng, n) -> array of n outcomes. Everything else - parallelizing across every core, keeping memory flat no matter how many trials run, tracking convergence, checkpointing a long run so it survives an interruption, and turning the result into a report worth looking at - is the engine's job.

from ulam import Simulator

def dart_throw(rng, n):
    x, y = rng.uniform(-1, 1, n), rng.uniform(-1, 1, n)
    return (x**2 + y**2 <= 1.0).astype(float) * 4.0

result = Simulator(dart_throw, seed=42).run(n_samples=1_000_000_000)
print(result.summary())

n              = 1,000,000,000
mean           = 3.14157
std            = 1.6422
std error      = 5.1931e-05
95% CI         = [3.14146, 3.14167]
elapsed        = 15.09s  (66,289,987 samples/s)

One billion trials, accurate to five significant figures of pi, in fifteen seconds on an 8-core laptop.

Why this is harder than it sounds

A Monte Carlo loop that runs a million times is a for-loop. One that's meant to run billions of times, on real problems, hits three walls a toy script doesn't:

1. Memory. A billion floats is 8 GB just to hold the raw samples - before you've computed anything from them. ulam never stores the trials; every batch is folded into running aggregates (mean, variance, a histogram, a random subsample for percentiles) via a numerically stable parallel Welford update, then discarded.
2. Randomness across cores. Split the work across 8 processes the naive way (seed each worker with its PID, or the time) and you get correlated random streams - a bug that silently corrupts results without ever raising an error. ulam uses NumPy's SeedSequence.spawn() mechanism so every batch, in any process, gets a provably independent stream, and a run with the same seed reproduces bit-for-bit regardless of how many workers touched it.
3. Time. A pure-Python per-sample loop tops out around a few hundred thousand iterations per second. ulam requires the model to be vectorized (model(rng, n) returns all n outcomes as a NumPy array in one call), which is what makes 60-90 million samples/second on ordinary hardware possible in the first place.

What's inside

Module	What it does
ulam.engine	Simulator - the orchestrator: batch planning, workers, convergence, checkpointing
ulam.streaming_stats	Mergeable online statistics: Welford moments, a fixed-bin histogram, a P² quantile estimator, a bounded random-subsample accumulator
ulam.rng	Reproducible parallel RNG streams (StreamManager) and an AntitheticGenerator for automatic variance reduction
ulam.sampling	Variance-reduction primitives: antithetic variates, Latin hypercube sampling, stratified sampling, control variates, importance-sampling weights
ulam.parallel	Batch scheduling and execution across processes/threads
ulam.convergence	Standard-error tracking and early stopping once a target precision is reached
ulam.checkpoint	Save/resume state for runs that take long enough to want that
ulam.report	Self-contained, offline-viewable interactive HTML reports (Plotly inlined)
ulam.cli	ulam run <model> --n <count> ... with a live terminal progress panel
ulam.models.*	16 ready-to-run example models (below)

Install

cd ulam
python3 -m venv .venv && source .venv/bin/activate
pip install -e ".[all]"     # numpy + scipy core, plus plotly/jinja2 (reports), rich (CLI), pytest (dev)

numpy and scipy are the only hard requirements; plotly/jinja2 (reports), rich (pretty CLI progress), and numba (optional JIT for user models) are opt-in extras so the core engine stays lightweight.

Command line

ulam list                                                  # see all 16 built-in models
ulam run pi --n 1_000_000_000                              # a billion trials, all cores, live progress
ulam run option_call --n 10_000_000 --report report.html   # write an interactive HTML report
ulam run craps --n 5_000_000 --seed 42 --workers 4         # reproducible, pinned worker count
ulam run --module my_model.py --func simulate --n 1e8      # bring your own model

The 16 built-in models

Every one of these is a real, checkable simulation - not a toy - and every one is picklable, so every one actually runs in parallel across all your cores (a surprisingly common failure mode for naive multiprocessing + closures, which ulam's model registry is built specifically to avoid).

Name	What it simulates
pi	Dart-throw estimate of pi (the canonical first Monte Carlo example)
buffon	Buffon's needle - an 18th-century geometric-probability estimate of pi
dice_sum	Sum of two six-sided dice
coin_streak	Longest heads streak in 100 coin flips
craps	A fully vectorized pass-line craps engine - validates against the known ~49.29% house edge
option_call / option_put	European options under geometric Brownian motion - validated against closed-form Black-Scholes
asian_option	Path-dependent Asian option (no closed form - has to be simulated)
portfolio_var	Correlated multi-asset portfolio risk (feeds into Value-at-Risk / Expected Shortfall)
random_walk / walk_2d	1D and 2D random walk endpoint distributions
first_passage	Steps for a 1D random walk to hit a barrier
project_risk	PERT/CPM critical-path project-completion risk, given a task network
epidemic_size / epidemic_peak	Stochastic chain-binomial SIR epidemic model - final size and peak load
queue_wait	M/M/1 queue mean wait time via the Lindley recursion - validated against queueing theory

Each is a few dozen lines in ulam/models/, meant to read as a template for writing your own.

Bring your own model

The entire contract is one function:

def my_model(rng: np.random.Generator, n: int) -> np.ndarray:
    """Return n independent outcomes, vectorized."""
    ...

from ulam import Simulator
result = Simulator(my_model, seed=42, n_workers=8).run(n_samples=500_000_000)

For multiprocessing, my_model needs to be a module-level function or a picklable callable object (a class with __call__) - the same rule as Python's multiprocessing in general. Every built-in model that needs internal state (a portfolio's covariance matrix, an epidemic's parameters) is written as a small @dataclass with __call__ for exactly this reason - a closure that captures those parameters would silently break the moment it's sent to a worker process. If your model genuinely can't be made picklable (an interactive notebook lambda, say), pass n_workers=1 and it still runs correctly, just single-threaded.

Variance reduction, for free or nearly free

Simulator(my_model, antithetic=True).run(n_samples=10_000_000)

Wraps the RNG your model receives so that calls to .random(), .uniform(), .standard_normal(), and .normal() automatically draw in mirrored antithetic pairs (u, 1-u) / (z, -z). For a model whose response is monotonic in its random inputs - which covers a large share of financial and physical models - this is close to free variance reduction. On the built-in European call option it delivers a measured 3.6x variance reduction (same accuracy for a third of the samples, or three times tighter for the same sample count) with a one-line change. ulam.sampling also provides Latin hypercube sampling, stratified sampling, and control variates for building custom variance reduction into your own model.

Convergence-aware runs

result = Simulator(my_model).run(
    n_samples=10_000_000_000,        # a hard ceiling...
    target_relative_se=0.001,        # ...but stop as soon as we're within 0.1% precision
)

Standard error shrinks like 1/sqrt(n) - quadrupling your sample count only halves your uncertainty. Past a certain point, most models have already told you what they're going to tell you, and running the rest of a ten-billion-trial ceiling is wasted compute. target_relative_se / target_absolute_se check convergence after every batch and stop the instant the target is hit, which is what actually makes an enormous n_samples a safe default to reach for rather than a knob you have to hand-tune.

Checkpointing

result = Simulator(my_model).run(
    n_samples=5_000_000_000,
    checkpoint_path="run.ulamckpt",
    checkpoint_every_s=30,
)

If the process is killed partway through, re-running the identical call resumes from the last checkpoint rather than starting over - because every batch's random stream is a pure function of (seed, batch index), a resumed run consumes exactly the random draws the original would have, batch for batch. A run that finishes cleanly deletes its own checkpoint file.

Reports

from ulam.report import render_report
render_report(result, "report.html", model_name="option_call")

Produces one self-contained HTML file - Plotly's JS is inlined, so it's safe to open offline or attach to an email - with an interactive histogram, a running-mean-with-confidence-band convergence chart, a throughput chart, and a percentile table, styled as a dark instrument panel rather than a wall of printed numbers. This is what ulam run <model> --report out.html writes.

See reports/showcase_asian_option.html for a live example - an Asian option priced with antithetic variance reduction, converging early once it hit its precision target.

Performance

Measured on an 8-core machine, the pi model (two uniform draws and a comparison - about as cheap as a model gets):

Trials	Workers	Wall time	Throughput
2,000,000	1	0.1s	~18M samples/s
500,000,000	8	~6-12s	~45-85M samples/s
1,000,000,000	8	~11-23s	~45-90M samples/s

(Throughput varies with system load; run python benchmarks/benchmark.py to measure your own machine.) Heavier models - path-dependent options, epidemic simulations with hundreds of internal time steps - cost proportionally more per trial, but the scaling story is identical: vectorized batches, parallelized across cores, streaming aggregation so memory never grows with n.

Tests

pytest tests/ -v      # 43 tests: streaming statistics, engine correctness, all 16 models, RNG/sampling, reports, CLI

Tests validate against ground truth wherever one exists - the option pricer against closed-form Black-Scholes, the craps engine against the known house edge, the M/M/1 queue against queueing theory's closed-form expected wait - not just "it runs without crashing."

Project layout

ulam/
  engine.py            Simulator, SimulationResult
  streaming_stats.py   RunningMoments, RunningHistogram, P2Quantile, Accumulator
  rng.py               StreamManager, AntitheticGenerator
  sampling.py           antithetic/LHS/stratified sampling, control variates
  parallel.py          batch scheduling, process/thread execution
  convergence.py       ConvergenceTracker, early stopping
  checkpoint.py        CheckpointState (save/resume)
  report.py            interactive HTML report generation
  cli.py               command-line interface
  models/              16 built-in example models
tests/                 43 tests across every module
examples/              quickstart, a full billion-trial run, a custom model walkthrough
benchmarks/            throughput benchmark script

License

MIT.