Lesson 09 - Engineering the Simulation¶

The beauty of Monte Carlo simulation is in its simple idea: repeat many simple, possibly random steps, and average the results. However, as a simulation grows in complexity, it naturally becomes a small piece of engineering software. Whether its users are many or one, good software requires documentation and some way to ensure reliable outputs. This lesson is about the habits that make Project 1 easier to build, debug, test, and trust.

By the end of this lesson, you should be able to:

separate exploratory notebook work from reusable Python functions;
design small functions with clear inputs, outputs, units, and assumptions;
write unit tests for deterministic pieces of a simulation;
use defensive programming to catch impossible inputs and runaway histories;
make random simulations reproducible without hiding uncertainty;
use a debugger to inspect program state instead of relying only on print statements;

Introduction¶

The previous lessons introduced histories, random variables, sampling, estimators, and uncertainty. Those ideas answer questions such as:

How do we sample a physical event?
How do we estimate a mean response?
How does the standard error change with the number of histories?
How do we interpret variation between independent runs?

Those questions are fundamentally mathematical, but can we implement the mathematics reliably?

For example, suppose a simulation estimates detector response from a slab. The statistical result might look reasonable, but the code could still be wrong. A reasonable-looking answer can come from a wrong distance-to-boundary check, a missed absorption branch, a collision counted in the wrong bin, or a temporary variable left over from yesterday’s debugging session. Monte Carlo methods are especially vulnerable to this because sampling noise can hide small mistakes.

Engineering the simulation means making the code answerable. You should be able to point to a function and say:

what it assumes;
what it returns;
how it handles invalid input;
how it was tested;
what random number generator it uses;
how you would debug it if it behaved strangely.

A useful simulation is not just one that runs. It is one that you can trust through verification.

The next block of the course moves toward deterministic balance equations and linear systems. Sometimes simulation gives an estimate by sampling histories. Sometimes the same physical thinking can be written as a balance equation and solved directly. In both cases, the discipline of careful structure absolutely matters.

The simulation-engineering workflow¶

A practical workflow for programming projects looks like this:

read the project statement
        ↓
complete one function at a time
        ↓
run given unit tests
        ↓
run your sanity checks
        ↓
check results and uncertainty
        ↓
communicate conclusions

Each stage has a different job.

This separation feels like extra work at first, but it saves time as soon as a project has more than a few functions.

A small Project 1 layout¶

A simple layout for project 1 is

path/to/project1/
  lastname_firstname_project1.py
  ne415_project1_student_test.py
  lastname_firstname_collision_profile.png
  optional_exploration.ipynb

Here, /path/to might be C:\Users\jaroberts\Documents\ne415\project1.

Keep the unit test file in the same folder as your completed project file while you work. The notebook is optional but encouraged, especially for communicating. The required simulation logic belongs in the .py file so that it can be imported and tested.

If your file is named doe_jane_project1.py, the test file should import that file. Then it can call functions such as:

sample_free_path(Sigma, rng)
sample_source(d, rng, source_type="boundary")
simulate_history(d, Sigma, c, f, bin_edges, rng, x0=0.0, direction0=+1)
mean_and_se(samples, axis=0)

That is the most important organizational rule for Project 1: the functions described in the starter file must exist, must have the expected names, and must return the expected types.

Keep notebooks tidy and reproducible¶

Notebooks are excellent for engineering analysis. They let you combine text, equations, code, plots, and interpretation in one document. They are not as good for storing the core logic of a simulation.

A notebook becomes fragile when it depends on hidden state: variables created ten cells above, a function redefined halfway through, a random seed set in a cell that was not rerun, or an array modified in place without being obvious. The longer the notebook becomes, the harder it is to know whether the visible output could be reproduced from a clean start.

Tidy habits¶

A good notebook for a simulation project usually does four things reproducibly:

defines or imports parameters;
calls functions from your completed project file;
makes plots and tables;
explains the results.

Before accepting notebook results, do a clean run:

Kernel → Restart Kernel and Run All Cells

Then check whether the same conclusions appear. This is one of the simplest reproducibility tests you can do.

Other good habits include:

keep parameter definitions near the top
move reusable code into your Project 1 .py file
restart and run all before saving final outputs
record the random seed used for any displayed result
avoid long outputs
again, restart and run all!

Useful Jupyter magics¶

Jupyter and IPython provide special commands called magics that provide tidy ways to enhance notebooks. A line magic starts with % and applies to one line. A cell magic starts with %% and applies to an entire cell.

The line magics include:

%pwd          # print the current working directory
%ls           # list files in the current directory
%time expr    # time one expression
%lsmagic      # list available magics

The cell magics include:

%%time
# time the whole cell

%%writefile path/to/file.py
# write the contents of this cell to a file

Design functions before writing long simulations¶

A simulation history is a sequence of small operations. Each operation should be as clear and testable as possible.

Consider a one-dimensional slab model. Even before writing the whole history loop, we might need functions for:

validating material parameters;
sampling a free path;
computing distance to the nearest boundary;
deciding whether a collision absorbs or scatters;
updating position after a flight;
recording the outcome of a history;
computing the mean and standard error of many histories.

Each of those pieces is easier to understand than the full simulation. Many of them can be tested with known inputs.

Give each function one job¶

A function should have a job that can be described in one sentence. For example:

sample_free_path(Sigma, rng)

Samples one exponentially distributed free path for a positive total macroscopic cross section.

sample_source(d, rng, source_type="boundary")

Returns an initial source position and direction. For the boundary source, this is \(x_0 = 0\) and direction \(+1\). For the volume source, this samples \(x_0\) uniformly in the slab and samples direction isotropically.

mean_and_se(samples, axis=0)

Computes the sample mean and estimated standard error of the mean for detector tallies or collision-bin tallies.

Each of those functions has clear inputs and outputs. Each can fail loudly when given impossible inputs. Each can be tested separately.

Prefer pure functions when possible¶

A pure function returns a value based only on its inputs and does not change outside state. Pure functions are easier to test because the same input always gives the same output.

This helper is pure:

[1]:

import numpy as np


def free_path_from_uniform(Sigma, xi):
    """Transform one uniform random number into an exponential free path."""
    if Sigma <= 0.0:
        raise ValueError("Sigma must be positive")
    if not (0.0 < xi < 1.0):
        raise ValueError("xi must satisfy 0 < xi < 1")

    return -np.log(xi) / Sigma

There is no randomness here. There is no plotting. There is no hidden global variable. That makes the function easy to test.

A function that samples from a random distribution is not pure because it consumes random numbers. In Project 1, the public function can stay very small:

[2]:

import numpy as np


def sample_free_path(Sigma, rng):
    """Sample a free path using a NumPy random Generator."""
    return free_path_from_uniform(Sigma, rng.random())

Now free_path_from_uniform can be checked exactly for selected values of xi, while sample_free_path is a tiny wrapper around the random number generator.

Make units and assumptions visible¶

Many engineering bugs are unit bugs. A function name such as sample_free_path is helpful, but the docstring should also say what units are expected. If \(d\) and position variables such as \(x\) and \(x_0\) are in centimeters and \(\Sigma\) is in inverse centimeters, then the sampled free path is in centimeters.

A good docstring does not need to be long. It should answer:

what the function computes;
what each input means;
what units are expected when units matter;
what the function returns;
what exceptions are raised for invalid input.

This is not just for other people. It is for you two weeks later.

Unit tests: small checks with known answers¶

A unit test checks one small piece of code. It should be fast, focused, and deterministic. The built-in Python module unittest is enough for this course.

A typical unit test follows the pattern arrange, act, assert:

arrange the inputs;
act by calling the function;
assert that the output has the expected behavior.

The provided student test file follows this same pattern. For example, a sampling test can force a known random number and check the expected free path:

import unittest
import numpy as np

import lastname_firstname_project1 as project1


class FakeRNG:
    def __init__(self, values):
        self.values = list(values)

    def random(self):
        return self.values.pop(0)


class TestProject1Sampling(unittest.TestCase):

    def test_free_path_known_random_number(self):
        rng = FakeRNG([np.exp(-1.0)])
        self.assertAlmostEqual(project1.sample_free_path(Sigma=2.0, rng=rng), 0.5)


if __name__ == "__main__":
    unittest.main()

After editing the import line in ne415_project1_student_test.py, run the tests from the folder containing your project file:

python ne415_project1_student_test.py

The test output is useful while learning because the test names become a checklist of expected behavior.

Test names should explain behavior¶

A test named test_1 is not helpful. A test named test_invalid_position_raises says exactly what should happen. When a test fails, the name should help you understand the intended behavior before you even open the file.

Good test names often follow this pattern:

test_<condition>_<expected_behavior>

Examples:

test_empty_input_raises
test_known_three_value_mean
test_zero_histories_raises
test_sample_source_volume_position
test_sample_source_volume_direction
test_sample_source_boundary_position
test_sample_source_boundary_direction
test_scatter_probability_above_one_raises

Test floating-point values carefully¶

Do not usually test floating-point results with exact equality. Roundoff is part of numerical computing. Use assertAlmostEqual for scalar values or np.allclose for arrays.

For example:

[3]:

import unittest
import numpy as np

class TestSampling(unittest.TestCase):

    def test_free_path_known_uniform(self):
        # If xi = exp(-1), then the exponential distance is 1 / Sigma.
        xi = np.exp(-1.0)
        self.assertAlmostEqual(free_path_from_uniform(Sigma=2.0, xi=xi), 0.5)

This test checks the transformation without depending on the random number generator.

Test exceptions, not just successful cases¶

Invalid inputs should fail loudly. That behavior should also be tested.

[4]:

import unittest
import numpy as np

class TestSamplingValidation(unittest.TestCase):

    def test_Sigma_must_be_positive(self):
        with self.assertRaises(ValueError):
            free_path_from_uniform(Sigma=0.0, xi=0.5)

    def test_xi_must_be_between_zero_and_one(self):
        with self.assertRaises(ValueError):
            free_path_from_uniform(Sigma=1.0, xi=1.0)

A test like this may seem unimportant until an invalid parameter silently contaminates a long simulation. Then it becomes one of the most important tests in the project.

Turn bugs into tests¶

When you find and fix a bug, write a test that would have caught it.

A good workflow is:

find bug → reproduce bug with a small input → write failing test → fix code → watch test pass

The failing test is valuable because it proves you understand the bug. The passing test is valuable because it protects future edits.

Defensive programming: fail early and clearly¶

Defensive programming means writing code that catches impossible or suspicious situations close to where they occur. It is better for a simulation to stop with a clear error message than to continue for ten minutes and produce a plausible but meaningless result.

Validate inputs at the boundary¶

Every public-facing function should check its most important assumptions. A function that runs histories might validate:

[5]:

def validate_simulation_inputs(N, d, Sigma, c, f, M, max_collisions):
    if N <= 0:
        raise ValueError("N must be positive")
    if d <= 0.0:
        raise ValueError("d must be positive")
    if Sigma <= 0.0:
        raise ValueError("Sigma must be positive")
    if not (0.0 <= c <= 1.0):
        raise ValueError("c must satisfy 0 <= c <= 1")
    if not (0.0 <= f <= 1.0):
        raise ValueError("f must satisfy 0 <= f <= 1")
    if M <= 0:
        raise ValueError("M must be positive")
    if max_collisions <= 0:
        raise ValueError("max_collisions must be positive")

The checks are not complicated. Their value is that they move the failure to the beginning of the run. If c = 1.2, you want to know immediately.

Use specific exceptions and helpful messages¶

This is not helpful:

raise Exception("bad")

This is better:

raise ValueError("c must satisfy 0 <= c <= 1")

The exception type tells you what kind of problem occurred. The message tells you what assumption was violated.

Some useful built-in exception types are:

Exception type	Use when
`ValueError`	An argument has the right type but an invalid value, such as `Sigma <= 0` or `c > 1`.
`TypeError`	An argument has the wrong type, such as passing a string where a number is required.
`IndexError`	Code tries to access an array or list index that does not exist.
`KeyError`	Code tries to access a missing dictionary key, such as a missing result field.
`RuntimeError`	The program reaches a state that should not happen during normal execution, such as exceeding `max_collisions`.
`FloatingPointError`	A numerical value becomes non-finite or invalid in a place where that should be treated as a failure.

Guard against runaway histories¶

Random walks can sometimes take many steps. A bug can make them take infinitely many steps. Use a maximum step count or collision count to protect the simulation.

A loop like this is dangerous:

while alive:
    # update particle state
    ...

If alive never becomes False, the history never ends. A safer pattern is:

for collision in range(max_collisions):
    # update particle state
    if history_has_ended:
        return result

raise RuntimeError("history exceeded max_collisions without terminating")

The for loop makes the limit visible. The RuntimeError says that the simulation reached a state that should not happen under normal assumptions.

Check for non-finite values¶

Numerical bugs often appear as nan or inf. These can spread through arrays silently.

if not np.isfinite(distance):
    raise FloatingPointError(f"non-finite distance: {distance}")

This check is especially useful inside a history loop while debugging. You may later remove or reduce expensive checks once the code is stable, but early in a project it is better to catch problems close to their source.

Assertions are for internal logic, not user input¶

An assert statement is useful for internal invariants during development:

assert collisions >= 0
assert outcome in {"absorbed", "leaked_left", "detected"}

Use explicit raise ValueError(...) for validating user inputs or project parameters. Assertions can be disabled when Python runs in optimized mode, so they should not be the only protection against invalid external input.

Random number generators and reproducibility¶

Modern NumPy code should create a random number generator object and pass it into functions that need randomness:

[6]:

import numpy as np

rng = np.random.default_rng(202609)
value = rng.random()

Then simulation functions can accept rng as an argument:

def simulate_history(d, Sigma, c, f, bin_edges, rng, x0=0.0, direction0=+1):
    ...

This is better than having functions call the global np.random state internally. Passing the generator makes randomness visible in the function interface. It also makes tests easier because you can pass a seeded generator or a fake generator.

Seed once near the top-level run¶

A good habit is to create one generator for a simulation run, then pass it downward. In the Project 1 starter, run_simulation accepts a seed, so run_simulation is the natural place to create the generator:

def run_simulation(N, d, Sigma, c, f, M, seed=None,
                   source_type="boundary"):
    rng = np.random.default_rng(seed)
    bin_edges = np.linspace(0.0, d, M + 1)

    for i in range(N):
        x0, direction0 = sample_source(d, rng, source_type)
        detector_tally, collision_tally = simulate_history(
            d, Sigma, c, f, bin_edges, rng, x0=x0, direction0=direction0
        )
        ...

Then a demonstration run passes the seed into run_simulation:

results = run_simulation(N=50_000, d=5.0, Sigma=1.0,
                         c=0.6, f=0.5, M=20, seed=202609)

Avoid reseeding inside low-level functions. This is a common mistake:

[7]:

def sample_free_path(Sigma):
    np.random.seed(12345)       # Do not do this inside the sampler.
    return np.random.exponential(1.0 / Sigma)

Reseeding inside a function can accidentally produce the same “random” value over and over. It also makes the function interfere with other code that uses randomness.

Record enough information to rerun the case¶

For every result you plan to report, record:

the model parameters;
the number of histories;
the random seed or seed sequence;
the estimator definition;
the date or version of the code if available;
the final estimate and uncertainty.

A simple result dictionary can help:

result = {
    "N": N,
    "d_cm": d,
    "Sigma_per_cm": Sigma,
    "c": c,
    "f": f,
    "M": M,
    "seed": seed,
    "detector_mean": results.detector_mean,
    "detector_se": results.detector_se,
}

A dictionary like this can be saved as JSON or converted into a table. It also makes it easier to label plots correctly.

Reproducible does not mean certain¶

Monte Carlo results have sampling uncertainty. If you run the same code with the same seed, you should get the same result. If you run the same code with a different seed, you should expect a slightly different result. That difference is not automatically a bug. It is part of the method.

A good project report should distinguish:

reproducibility: can I rerun this exact calculation?
uncertainty: how much random variation should I expect?
correctness: did I implement the intended model?

Those are related, but they are not the same.

A practical debugging workflow¶

When a simulation behaves strangely:

Save the strange case. Record the parameters and seed.
Make it small. Try fewer histories or a simplified input that still reproduces the problem.
Move the relevant code into your Project 1 file. Debugging a .py file is usually cleaner than debugging a long notebook.
Run one tiny case. Call one function with a small, controlled input.
Set breakpoints. Pause at the start of the function or before the suspicious branch.
Step through. Watch x, mu, sampled distances, event choices, counters, and outcomes.
Add a test. Once you understand the bug, write a regression test.

A tiny debug driver might be:

import numpy as np

rng = np.random.default_rng(12345)
bin_edges = np.linspace(0.0, 5.0, 11)

history = simulate_history(
    d=5.0,
    Sigma=1.0,
    c=0.6,
    f=0.5,
    bin_edges=bin_edges,
    rng=rng,
    x0=0.0,
    direction0=+1,
)

print(history)

Open the project folder in VS Code, set breakpoints inside simulate_history, and run your file or a small scratch cell under the debugger.

What to inspect¶

For a particle-history simulation, useful variables to watch include:

current position;
current direction;
sampled distance to collision;
distance to boundary;
scatter and direction-change probabilities;
random values used for event choices;
collision count;
detector score and collision-bin tallies;
termination condition.

The goal is to compare what the code does with what the algorithm says it should do. When those differ, you have found the neighborhood of the bug.

Conditional breakpoints¶

A conditional breakpoint pauses only when a condition is true. This is useful when a bug appears after many iterations.

Examples:

collision == 17
x < 0
not np.isfinite(track_length)
outcome == "absorbed"

Instead of printing thousands of lines, let the debugger stop when the suspicious condition occurs.

Read tracebacks from bottom to top¶

When Python raises an exception, the traceback shows the call stack. The bottom of the traceback usually shows the exact line where the error occurred. The lines above show how the program got there.

Do not ignore tracebacks. They are maps. Start near the bottom, identify the line that failed, then move upward to understand which function called it with which inputs.

Restart the kernel when debugging notebooks¶

If a notebook “makes no sense,” suspect hidden state. Restart the kernel and rerun the notebook. Also check that the notebook imported the project file you think it imported:

import lastname_firstname_project1
print(lastname_firstname_project1.__file__)

Many mysterious notebook bugs are actually stale definitions or path problems.

Flowcharts: making the algorithm visible¶

A flowchart is a picture of control flow: sequence, branches, loops, and termination. It is not just decoration. It helps you see whether the algorithm has all necessary exits.

Common flowchart shapes are:

Start/end:      ([text])
Process:        [text]
Decision:       {text?}
Input/output:   [/text/]

Mermaid is a lightweight text format for diagrams. A simple loop looks like:

flowchart TD
  S([Start]) --> I[i = 0]
  I --> C{i < n?}
  C -- yes --> W[do work]
  W --> U[i = i + 1]
  U --> C
  C -- no --> E([End])

A particle-history flowchart might look like this:

flowchart TD
  S([Start history]) --> P[/Read parameters and RNG/]
  P --> V{Inputs valid?}
  V -- no --> VE[Raise ValueError]
  V -- yes --> I[Initialize x, direction, tallies]
  I --> C{collision count < max?}
  C -- no --> RE[Raise RuntimeError]
  C -- yes --> F[Sample free path]
  F --> B[Compute distance to boundary]
  B --> L{Boundary reached first?}
  L -- detector --> O[Return detected result]
  L -- left leak --> LL[Return leaked-left result]
  L -- no --> M[Move to collision point]
  M --> A{Absorbed?}
  A -- yes --> R[Return absorption result]
  A -- no --> D[Sample scatter direction]
  D --> C

Use the flowchart to ask practical questions:

Does every loop have a way to stop?
Are invalid inputs caught before the loop begins?
Are all possible outcomes returned?
Are there branches for detection and left leakage before collision?
Is there a maximum-collision guard?
Is scoring separate from transport logic?

Before writing a complicated history loop, sketch the flowchart. After writing the loop, compare the code to the chart. If they differ, decide which one is wrong.

Repeatable runs¶

A notebook is good for exploration and communication. The main() function in your Project 1 file is good for a repeatable demonstration run. The important guard is:

if __name__ == "__main__":
    main()

That guard lets the student test file import your functions without accidentally running a long simulation. Inside main(), keep the case clear and reproducible:

def main():
    seed = 202609
    N = 50_000
    d = 5.0
    Sigma = 1.0
    c = 0.6
    f = 0.5
    M = 20

    results = run_simulation(
        N=N, d=d, Sigma=Sigma, c=c, f=f, M=M, seed=seed,
        source_type="boundary",
    )

    print(f"Detector response = {results.detector_mean:.6f} +/- {results.detector_se:.6f}")
    plot_collision_profile(results, d, Sigma, c, f, N, "collision_profile.png")

if __name__ == "__main__":
    main()

This gives you one command that reruns a known case:

python lastname_firstname_project1.py

Documentation habits for student projects¶

Documentation does not need to be fancy. It needs to preserve decisions.

At minimum, a project should explain:

what physical or mathematical model is being simulated;
what assumptions are made;
what each main parameter means;
what units are used;
what estimator is reported;
how uncertainty is estimated;
how to run the demonstration case;
how to run the tests;
what files contain the core code.

A short comment block near main() can be enough:

# Demonstration case:
# - boundary source at x = 0 cm with initial direction +1
# - d = 5 cm, Sigma = 1 cm^-1, c = 0.6, f = 0.5
# - N = 50,000 histories, M = 20 collision bins
# - outputs detector response and collision_profile.png

Notes like this make the project easier to grade, easier to revisit, and easier to debug.

Comments should explain why¶

Do not comment every line. Comments are most useful when they explain why something is done, not what Python syntax already says.

Less useful:

x = x + mu * s  # add mu times s to x

More useful:

# Move only to the collision site; boundary crossings were already handled above.
x = x + mu * s

The first comment repeats the code. The second comment explains the algorithmic assumption.

Bridge to Lesson 10¶

Monte Carlo simulation estimates behavior by sampling histories. That is powerful because it can handle complicated random processes. It also means that code structure, testing, reproducibility, and uncertainty matter.

In the next lesson, we ask a different question:

Sometimes simulation gives us an estimate by sampling histories. What changes when the model can be written as a balance equation and solved directly?

That question leads to linear systems of the form

\[A x = b.\]

The programming habits from this lesson still apply. We will still separate exploratory notebooks from reusable functions. We will still check residuals and validate inputs. We will still prefer clear structure over one giant block of code. The mathematical tool changes, but the engineering habit remains the same: build models that can be inspected, tested, and trusted.

References and further reading¶

These are optional references for the tools used in this lesson:

Python documentation: `unittest — Unit testing framework <https://docs.python.org/3/library/unittest.html>`__
NumPy documentation: Random ``Generator` <https://numpy.org/doc/stable/reference/random/generator.html>`__
Scientific Python Blog: Best Practices for Using NumPy’s Random Number Generators
IPython documentation: Built-in magic commands
Visual Studio Code documentation: Python debugging in VS Code