Skip to main content

Parallel Simulation Experiments on a Personal Computer

7 Practical Considerations

The earlier chapters gave you a working formula. You write a worker function, hand it to sim_template, pick a number of cores, and collect the results. For a great deal of real work that is genuinely all you need. This chapter is about the gap between code that works once on your machine and code you can trust in the long run, hand to a collaborator, or come back to in six months and still understand.

None of what follows is deep; but, most of it is the kind of thing you will only discover after the fact:

  • How many cores do you actually have, as opposed to how many the operating system claims?
  • Why does adding cores sometimes make a simulation slower instead of faster?
  • What do you do when you stop a run and your laptop’s fans keep roaring an hour later?
  • How do you time code in a way that means something, name things so the next person can follow them, and seed a study so that “reproducible” is a fact rather than a hope?

The sections that follow tackle these questions one at a time. Section 7.1 is about counting cores honestly. Section 7.2 covers the bottleneck that is not speed at all but memory. Section 7.3 is about recovery when a parallel run goes wrong, including the orphaned processes that can outlive the program that spawned them. Section 7.4 replaces eyeballing a clock with proper benchmarking. Section 7.5 is about naming, comments, and project layout, the unglamorous habits that make parallel code survivable. Section 7.6 looks at load balancing, why the work does not always divide as evenly as the arithmetic suggests. And Section 7.7 closes with reproducibility, which is the whole point of fixing a seed in the first place. You can peruse this now to get an overview and come back to it when you need it.

7.1 Cores, Threads, and How to Count Them

Every parallel call in this book asks you for the same number, n_cores. The quality of that choice does more for your wall-clock time than almost anything else. To choose well you have to know what you are selecting from, and in this instance the machine can be less than forthcoming.

The obvious move is to ask the computer how many cores it has. In R that is

library(parallel)
detectCores()

and in Python it is

import os
os.cpu_count()

The trouble is that the number these commands return is not the number of cores you want to use, for two separate reasons. The first is hyperthreading. Many Intel and AMD chips present each physical core to the operating system as two logical processors. A genuine four-core chip reports eight. Hyperthreads share the underlying arithmetic hardware. Consequently, for the compute-bound numerical work in this book they deliver far less than what would be obtained from the addition of another physical core. Counting them as cores and setting n_cores to the reported number that includes hyperthreads is a common way to be disappointed.

The second reason is the one introduced back in Chapter 1: modern machines mix performance and efficiency cores. The two iMac and ASUS machines used for timing in Chapters 3 and 4 make the point. The iMac is an eight-core Apple M4 with four performance and four efficiency cores; the ASUS is a twenty-core Windows laptop with eight performance and twelve efficiency cores. Efficiency cores exist to save power on light tasks, not to crunch numbers, and detectCores cheerfully counts them alongside the rest.

There is also a subtler point that Chapter 1 raised that bears repeating here, because it changes how you read every one of these numbers. When you call makeCluster(n_cores) in R or open a Pool(n_cores) in Python, you are not seizing hardware cores directly. You are asking the operating system to start n_cores worker processes, and the OS then time-shares those processes across whatever physical cores it sees fit, alongside your browser, your editor, and everything else running. The mapping from process to core is the operating system’s to make, not yours. This is why asking for more workers than you have performance cores rarely helps and often hurts: you have simply handed the scheduler more processes to juggle on the same hardware.

So how should you choose n_cores in practice? A reliable rule of thumb is to count the performance cores, not the logical-processor total, and start there, leaving one core free for the operating system and your own interactive work. On the iMac that means around four; on the ASUS, around eight. You do not have to guess these values. R’s parallelly package, a companion to the popular future framework, exposes parallelly::availableCores(), which is more conservative than detectCores() and respects environment limits set by, for example, a cluster scheduler. But the most honest answer comes from the machine itself: run the short speed-up experiment of Section 7.4, watch where the curve flattens, and let that be your n_cores. The number the hardware reports is a starting hypothesis; the number that stops getting faster is the truth.

7.2 Memory: The Other Bottleneck

It is natural to assume that the constraint on a simulation is always speed, and that more cores is therefore always the lever to pull. For a whole class of problems that assumption is wrong, and pulling the lever makes things worse. The other bottleneck is memory.

The reason is built into how the parallel tools in this book work. R’s default PSOCK cluster, and Python’s default process-spawning behavior on Windows and macOS, both start each worker as a fresh, independent process. A fresh process does not share your data; it gets its own copy. If your simulation reads a one-gigabyte dataset and you run it on eight workers, you can be holding eight gigabytes of it in RAM at once, plus the original, plus the operating system’s own needs. The arithmetic is unforgiving, and it is easy to walk into. A reviewer of an early draft of this book raised exactly this case: a residual bootstrap, of the kind we built in Chapter 5, runs not on a teaching dataset but on a ten-gigabyte production one. Two or three workers and you are out of memory before the processors have broken a sweat.

The symptom is distinctive once you have seen it, and it is the mirror image of the compute-bound case. When you are short on cores, adding cores makes things faster. When you are short on memory, adding cores makes things slower. Memory use climbs as you add workers, the machine begins swapping to disk to fake the RAM it does not have, throughput collapses because disk is thousands of times slower than memory, and in the worst case the operating system kills your run outright to save itself. If you ever see a study get slower as you grant it more cores, suspect memory before anything else.

A little estimation goes a long way here. Before launching a large run, find out how big the per-worker footprint actually is. In R, object.size on the data you are about to ship to each worker gives you the base figure; multiply by n_cores and compare against the RAM you have. In Python, sys.getsizeof is a rough first look and the memory_profiler package a more careful one. While a run is going, keep a system monitor open (Activity Monitor on macOS, Task Manager on Windows, top or htop on Linux) and watch the memory line rather than the CPU line. If it is the memory line that is pinned, more cores is the wrong fix.

The right fixes are the opposite of more parallelism. The simplest is to use fewer workers: three workers that fit in memory will finish a memory-bound job that eight workers will thrash to death. Beyond that, you can avoid the copying altogether. On macOS and Linux, R can build a FORK cluster with makeForkCluster() (and mclapply() uses the same mechanism), and Python’s multiprocessing can use the fork start method. A forked worker inherits the parent’s memory through the operating system’s copy-on-write machinery, which means a large dataset the workers only read is shared rather than duplicated, and stays shared until someone writes to it. That can turn an impossible eight-copy job into a comfortable one. The catch, and the reason this book defaults to PSOCK and spawn despite the cost, is portability: forking is not available on Windows and interacts badly with some multithreaded libraries, so it is a deliberate choice for a specific machine rather than a default to reach for blindly. When even forking is not enough, the honest conclusion is that the data has outgrown the laptop, which is the memory-bound signal that Chapter 8 addresses under when to move beyond a personal machine.

7.3 When Things Go Wrong: Debugging Parallel Code

Parallel code fails in ways serial code does not, and the failures are harder to see because the action is happening inside worker processes you are not watching. Two habits make the difference between a frustrating afternoon and a quick fix.

The first is to do all of your debugging in serial. Every sim_template in this book takes n_cores as an argument, and the single most useful thing you can do when something breaks is to set n_cores = 1. With one core, the worker runs in a context you can actually inspect: error messages come straight back to you instead of being swallowed and possibly re-reported by the cluster machinery, you can drop a print or a browser into the worker and see its output. The whole confusing layer of inter-process communication is simply absent. Get the worker correct and fast on a small problem with n_cores = 1, confirm it gives the answer you expect, and only then scale up the cores and the replications. A bug that survives to the parallel run is then almost always in the parallel plumbing, which is a much smaller place to look.

The second habit is to understand what happens when one worker experiences errors partway through a run. The frameworks differ in their responses. Python’s Pool.starmap will propagate an exception from any worker back to your main program, which stops the run. So a single bad replication takes the whole job down and you see the traceback. R’s parallel functions are quieter; depending on how you call them, an error in one task may surface as a cryptic message buried in the returned list rather than a clean stop. Either way the lesson is the same: validate the inputs and guard the edges of your worker before you launch ten replications. The cost of a crash on replication 9,998 is the whole run.

Then there is the failure the reviewers were most emphatic about, because it is the one that quietly costs you money and machine time: orphaned processes. When you build a cluster you are starting real, independent processes, and they do not automatically die when your script does. If you interrupt a run, close R or Python without a clean shutdown, or the controlling process crashes, those workers can keep running, holding CPU and memory, for hours or days, until you notice the computer’s fans whirring or the next run mysteriously has no cores to use. The discipline that prevents this is to make shutdown automatic rather than something you remember to do. In R, pair every makeCluster() with a stopCluster(), and put it where it will run even on failure:

cl <- makeCluster(n_cores)
on.exit(stopCluster(cl))   # runs even if the body errors out
# ... do the parallel work ...

The sim_template listings in Chapters 3 and 4 shut their clusters and pools down for exactly this reason; Python’s with Pool(...) as pool: form is the same idea, guaranteeing the pool is closed when the block exits no matter how that happens. When an orphan does still manage to get loose, you find it with the same system monitor from the previous section. Look for processes named R or python chewing CPU with no terminal attached, and end it from there: Activity Monitor or Task Manager has a button for it, and on macOS or Linux pkill -f sim_template or a plain kill on the process id will do it from the command line. Checking for strays after an interrupted run is a five-second habit that saves a great deal of confusion.

7.4 Benchmarking Properly

The simplest timing tool in R is system.time(). It’s the honest tool for the job: built in, needing no explanation, and for a run that takes seconds, good enough. For deciding things, though, a single system.time() reading is shaky ground. Run it twice and you will get two different numbers, because a real machine is never doing only your work: a background update kicks in, or the clock speed drifts under thermal limits. One stopwatch reading cannot tell signal from noise.

Proper benchmarking tools exist to resolve precisely this issue. In R the two standard choices are the bench and microbenchmark packages and we have used bench to obtain the median times, etc., that were reported for our R code. Both libraries run your expression many times and report a distribution, a median, and a spread, rather than a single number. They handle the details a hand-rolled timer gets wrong, such as warming up so the first slow run does not skew the result and accounting for garbage collection. A typical use looks like

library(bench)
bench::mark(
  serial   = sim_template(seed, n, n_reps, n_cores = 1),
  parallel = sim_template(seed, n, n_reps, n_cores = 4),
  # the two return equivalent results, not identical objects
  check = FALSE
)

Python’s equivalent ships with the language, the timeit module. This likewise runs a snippet repeatedly and reports the best of several timings, deliberately filtering out the slow outliers that come from the system doing something else. For longer simulation runs, where you care about wall-clock time rather than the microsecond precision timeit is built for, you can wrap time.perf_counter() around the call and that is perfectly adequate as long as you average a few runs.

What you do with these tools matters as much as the tools. The quantity worth measuring is speed-up, defined in Section 3.2 as the ratio of the serial time to the time on n_cores cores. Measure it not at one core count but across a range, one, two, four, eight, and look at the shape of the curve. Early on it climbs, perhaps not at the ideal rate but climbing. Then it bends over and flattens, and eventually, as the iMac showed in Section 3.2, it can turn back down. The core count where the curve stops climbing is your answer: it is the most parallelism this problem and this machine can usefully absorb, and paying for more, in cores, in cloud dollars, or in complexity, buys nothing. Benchmarking is not bookkeeping you do after the fact; it is how you find the right n_cores in the first place, and it costs a few minutes against runs that may take hours.

7.5 Naming, Code Style, and Project Organization

Simulation code has a way of outliving the afternoon it was written. The script that produced a figure gets reopened when a referee asks a question, handed to a student, or adapted for the next project, and at that point its readability stops being a matter of taste and becomes a matter of whether the work can be trusted and reused. A few habits here pay for themselves many times over.

Start with names. A reviewer of this book pushed back on calling everything simple_, on the grounds that “simple is in the eye of the beholder,” and the point generalizes. A name should describe an object’s purpose. We kept sim_template as our running example precisely because it describes its job, a template you fill in, rather than passing judgment on how hard that job is. We have tried to consistently name our worker functions as worker because they perform the same basic task: generate data and compute estimators. When you adapt this code, resist the temptation to call your worker f or your result x. The five seconds you save typing are borrowed against the hour someone later spends working out what x might contain.

Comment the seams, not the obvious. Ordinary serial code often needs little commentary; but, parallel code has joints that are genuinely not self-explanatory, and those are where a sentence earns its place: a note on why each worker gets its own random-number stream and what would break if they don’t, a reminder that the worker function must be self-contained because it runs in a process that has never seen your global environment, a line explaining how the replications are divided across cores. A collaborator reading your code can see what the parallel call does; the comment should tell them why it is safe.

Keep workers pure. The single most useful structural rule for simulation code is that the worker function should take its inputs, compute its result, and return it, with no side effects, no writing to files, no changing global state, no installing packages. An early version of this book’s own code installed a package from inside the worker function, and a reviewer rightly flagged it: a function that quietly alters your system while it runs is a function you cannot reason about, reproduce, or run safely on a shared machine. Installation belongs in your setup notes, run once, as the Preface describes; the worker should do only one thing.

Finally, give a study a shape on disk. A simulation project is naturally made of a few distinct kinds of file, and keeping them separate keeps the project legible: the worker function that defines the statistics, the sim_template driver that runs it in parallel, an analysis script that turns raw output into tables and figures, and a place for the results themselves. A layout as plain as

project/
  R/          worker functions and the sim_template driver
  analysis/   scripts that summarize and plot
  results/    saved output, named by the run that produced it
  data/       inputs, treated as read-only

is enough. The value is not tidiness for its own sake; it is that anyone, including you in six months, can see at a glance where the simulation is defined, how it was run, and where its output went. That is the difference between code that produced a result once and code that can produce it again.

7.6 Load Balancing

When sim_template divides n_reps replications across n_cores workers, it splits them as evenly as the arithmetic allows, giving each core a base share and parceling out any remainder one at a time, as the program-details sections of Chapters 3 and 4 described. If every replication took the same time, that would be the end of the matter and all cores would finish together. Real work is rarely so tidy, and the gap between an even division of tasks and an even division of time is what load balancing is about.

Two things spoil the tidy picture. The first is that the replications themselves may not be equal. In the segmented-regression study of Chapter 5, each replication runs a random search and an optimization whose running time depends on the data it happens to draw; some replications converge quickly and others grind. Hand a process a job that happens to contain several slow tasks and it will still be working long after its neighbors have gone idle, and the whole run is only as fast as that unluckiest process. The second is the hardware we have returned to throughout this chapter: on a machine that mixes performance and efficiency cores, an equal share of work handed to an efficiency core takes longer than the same share on a performance core, so equal division in this case guarantees uneven finishing times.

The lever that helps is chunk size, the trade between overhead and balance. Splitting the work into a few big chunks, one per process, keeps overhead low. This is because each worker is set up and torn down only once. But it balances badly, since one heavy chunk delays everything. Splitting into many small chunks, more tasks than processes, balances well, because a process that finishes early simply picks up the next waiting piece. But, each handoff carries communication cost, and too many tiny tasks drown the real work in dispatch overhead. There are tools that expose this directly. R’s parallel package offers load-balancing variants such as parLapplyLB and clusterApplyLB that hand out work dynamically rather than in one fixed split, and Python’s Pool accepts a chunksize argument and offers imap for the same purpose. The practical advice is unfussy: for the homogeneous problems that make up most of this book, the even split sim_template already performs is the right default and you need not think about it. When you notice, in the system monitor, that some cores go idle while others labor on, that is the signal to make the chunks smaller and let the fast cores help with the slow work.

7.7 Reproducibility

Every chapter of this book has handed sim_template a seed, and it is worth ending the practical material by being precise about what that seed does and does not buy you. Because reproducibility is the entire scientific justification for the care we have taken with random numbers. The promise from Chapter 2 is real: the same seed drives the same deterministic stream out of the generator, so a simulation run with a fixed seed gives bit-for-bit identical results every time it is run. That is what makes a simulated result a result at all, something a colleague can check rather than take on faith, and it is the reason the seed is the first argument to every template in the book rather than an afterthought.

There is one wrinkle that parallelism introduces and that you should know about before it surprises you. Because sim_template builds a separate, independent random-number stream for each process, as Chapter 2 and the program details of Chapter 3 lay out, the exact set of random draws depends on how many cores the work was split across. Run the same study with the same seed on four processes and then on eight, and you will get answers that are each perfectly reproducible but not identical to each other, because the second run divided the replications differently and therefore drew from the streams in a different pattern. This is not a bug and it does not threaten your conclusions, which are statements about a distribution that both runs estimate consistently. But it does mean the honest reproducibility recipe records n_cores alongside the seed. The full quartet, seed, n, n_reps, and n_cores, reproduces a run exactly; the seed alone reproduces it only on the same number of cores.

Reproducibility also extends past your own code to the software underneath it. A generator’s output, a package’s default, even a numerical routine can change between versions, so a study that ran identically last year may drift this year through no fault of your own. The defenses are simple and worth the small effort. Record the environment that produced a result: sessionInfo() in R and pip freeze or conda list in Python capture the package versions in a form you can save next to the output. Save the seed, the core count, and the parameters in or beside the results file itself, so a figure always carries the information needed to regenerate it. And where exact long-term reproducibility matters, pin your package versions rather than letting them float. A run you cannot reproduce is, scientifically, a run that did not happen.

7.8 Chapter Summary

In this chapter we stepped away from new simulations and toward the operational knowledge that makes the ones you have reliable. We saw that the core count a machine reports is a starting hypothesis rather than an answer, inflated by hyperthreads and efficiency cores, and that the performance-core count, confirmed by a speed-up experiment, is the number to trust. We met memory as the bottleneck that is not speed, the one where adding cores makes things slower, and the copying behavior of PSOCK and spawned workers that causes it, along with fewer workers and copy-on-write forking as the fixes. We covered recovery when a parallel run goes wrong: debugging in serial with n_cores = 1, knowing how each framework reports a worker’s error, and the orphaned processes that outlive an interrupted run unless you make shutdown automatic. We replaced a single stopwatch reading with proper benchmarking through bench, microbenchmark, and timeit, and with reading the speed-up curve to find where more cores stop helping. We argued for names that describe, comments on the parallel seams, pure side-effect-free workers, and a project laid out so it can be understood and rerun. We saw why an even split of tasks is not always an even split of time, and how chunk size trades overhead against balance. And we ended on reproducibility, the payoff for all the care with seeds, including the wrinkle that a parallel run reproduces exactly only when the core count is recorded alongside the seed.

These are the habits that separate code that worked once from code you can stand behind. With them in place, the natural next question is what to do when even a well-tuned laptop is no longer enough. That is where Chapter 8 turns next, to the cloud, the cluster, and the GPU.