Keyboard shortcuts

Press ← or → to navigate between chapters

Press S or / to search in the book

Press ? to show this help

Press Esc to hide this help

The Python package

The simulator ships as a Python package called tepsim. It is not a reimplementation and not a wrapper around a subprocess: it is the same Rust core the rest of this book documents, compiled into an extension module with PyO3. There is no C in the build, no Fortran, and no runtime dependency except NumPy.

A run started from Python is bit-identical to the same run started from Rust or from the browser, because all three call the same code with the same scenario. That is what makes the worked examples in the notebooks usable as evidence rather than as illustrations.

Installing it

It is not on PyPI. The version is still 0.0.0 and no release tag has been pushed, so nothing has ever been published. The publish job in .github/workflows/wheels.yml is wired up and gated on a v* tag, and the name is free, so pip install tepsim will be the line once there is a release. It is not the line today.

Install it from the repository instead:

$ pip install "git+https://github.com/jkitchin/tep-rust#subdirectory=crates/tepsim-py"

Two things about that command, both of which will bite otherwise.

It builds from source, so you need a Rust toolchain. There is no published wheel to fall back on, so pip clones the repository, invokes maturin, and compiles the whole core in release mode. That takes a couple of minutes on a warm machine and it fails partway through if cargo is not on your PATH. Install Rust from rustup.rs first. NumPy is pulled in automatically as a declared dependency, so nothing else is needed.

The quotes are not optional. Most shells treat # as the start of a comment, and without the quotes the #subdirectory= fragment is silently dropped, at which point pip tries to build the workspace root and fails with an error that does not mention the fragment at all.

The result is a cp39-abi3 wheel, so one build covers every GIL-enabled CPython from 3.9 upwards. Free-threaded interpreters have no stable ABI to target and need a version-specific build; 3.14t is the earliest that works, because PyO3 0.29 does not support the free-threaded build of anything below 3.14.

From a checkout

If you have the repository already, cargo xtask python is the supported path. It builds a release wheel, creates a throwaway virtualenv at .xtask-python/venv, installs the wheel and NumPy into it, proves that import tepsim resolves inside that virtualenv rather than somewhere else on the machine, and runs the binding's pytest suite against it.

$ cargo xtask python
$ .xtask-python/venv/bin/python -c "import tepsim; print(tepsim.__version__)"

That interpreter is the one the book's own tests use to check the quickstart transcript below. Note that cargo xtask python deletes and rebuilds the virtualenv every time it runs, so anything else installed into it, jupyter and matplotlib for the notebooks in particular, has to be reinstalled afterwards.

To build a wheel without installing it, maturin build --release -m crates/tepsim-py/Cargo.toml writes one and prints the path.

Quickstart

from concurrent.futures import ThreadPoolExecutor

import numpy as np

import tepsim as tep

print("tepsim %s: XMEAS(1..%d), XMV(1..%d), %d channels, IDV(1..%d)"
      % (tep.__version__, tep.MEASUREMENTS, tep.MANIPULATED, tep.CHANNELS,
         tep.DISTURBANCES))

# A Scenario says what to simulate, a Simulation does it, a Run holds what
# came out. That is the whole API.
run = tep.Simulation(tep.Scenario.baseline(seed=42, hours=48)).run()
print()
print("run:      %r" % run)
print("matrix:   %s %s" % (run.to_numpy().shape, run.to_numpy().dtype))
print("outcome:  %s" % run.outcome)
print("XMEAS(7): mean %.2f kPa over %.0f h"
      % (run.measurement(7).mean(), run.hours[-1]))

# The twenty disturbances say what they do, not only what the original header
# called them. Five of the published descriptions are the word "Unknown".
print()
print("the five the original leaves unexplained")
for fault in tep.faults():
    if fault.published == "Unknown":
        print("  IDV(%2d) %-9s %s" % (fault.index, fault.shape, fault.effect))

# Ground truth travels with the data, which the original records nowhere.
faulted = tep.Simulation(tep.Scenario.fault(1, hours=8)).run()
labels = faulted.labels()
print()
print("IDV(1) over 8 h: active at the last sample %s, %.2f h since onset"
      % (labels["active"][-1, 0], labels["since_onset"][-1, 0]))

# run() releases the GIL, so an ensemble is a thread pool and nothing has to be
# pickled.
sims = [tep.Simulation(tep.Scenario.fault(n, hours=8)) for n in range(1, 21)]
with ThreadPoolExecutor() as pool:
    runs = list(pool.map(tep.Simulation.run, sims))
print()
print("twenty 8-hour faulted runs: %d completed, %d tripped"
      % (sum(r.outcome == "completed" for r in runs),
         sum(r.outcome == "tripped" for r in runs)))

# A run is a pure function of its scenario, and a scenario is one line of text
# that parses back to an equal scenario.
scenario = tep.Scenario.fault(4, hours=8)
print()
print("digest:   %s" % scenario.digest)
print("text:     %s" % scenario.to_text())
print("parses back equal:      %s"
      % (tep.Scenario.from_text(scenario.to_text()) == scenario))
print("two runs bit-identical: %s"
      % np.array_equal(tep.Simulation(scenario).run().to_numpy(),
                       tep.Simulation(scenario).run().to_numpy()))
tepsim 0.0.0: XMEAS(1..41), XMV(1..12), 53 channels, IDV(1..20)

run:      <tepsim.Run 960 samples x 53 channels, 48.0 h, completed>
matrix:   (960, 53) float64
outcome:  completed
XMEAS(7): mean 2706.16 kPa over 48 h

the five the original leaves unexplained
  IDV(16) random    enables walk channel 9, the stripper steam valve capacity
  IDV(17) random    enables spike channel 10, the reactor coolant duty
  IDV(18) random    enables spike channel 11, the condenser coolant duty
  IDV(19) sticking  sticks valves 5, 7, 8 and 9; touches no equation in the model
  IDV(20) random    enables spike channel 12, the reactor outlet flow

IDV(1) over 8 h: active at the last sample True, 8.00 h since onset

twenty 8-hour faulted runs: 19 completed, 1 tripped

digest:   7f9accc04bb61e7f
text:     tepsim.scenario.v1;seed=4651207995;hours=8;step=2.777777777777778e-4;every=180;faults=4;controlled=1;idv12=0;trip=1;continuous=0;integrator=euler;events=
parses back equal:      True
two runs bit-identical: True

The API

Three classes carry the whole surface, and they divide the way the problem does. A Scenario is a description and holds no state. A Simulation is the machinery. A Run is the result.

Module level

NameValueMeaning
MEASUREMENTS41XMEAS(1..41)
MANIPULATED12XMV(1..12)
CHANNELS53a row: measurements then manipulated variables
DISTURBANCES20IDV(1..20)
DEFAULT_SEED4651207995.0the generator word compiled into teprob.f:1187
DEFAULT_STEP_HOURS1/3600one simulated second, the step INTGTR uses
DEFAULT_SAMPLE_EVERY180the three-minute spacing of d00 through d21
FORCED_DISTURBANCE_STEP28800the step the driver forces IDV(12) on at

channel_names() returns the 53 names in row order, so a CSV header and a matrix column cannot disagree about which is which. faults() returns the twenty Fault records.

Scenario

Immutable and cheap to copy, so a caller can build one, derive variants from it, and keep all of them.

There are four constructors: the bare Scenario(...), Scenario.baseline(...), Scenario.fault(n, ...) and Scenario.from_text(text). All but the last take the same keyword arguments, which are seed, hours, step_hours, sample_every, controlled, driver_forces_idv12 and trip_ends_the_run, plus faults on the bare one.

MemberWhat it is
seed, hours, step_hours, sample_everythe four numbers a run is measured in
controlledclosed loop under the published control scheme, or open loop with the valves held
driver_forces_idv12whether the driver switches IDV(12) on at hour eight, delta D-011
trip_ends_the_runwhether a shutdown stops the run, delta D-007
faultsthe active disturbances, one-based and ascending
steps, sampleshow many integrator steps and how many recorded rows
digesta content hash over everything affecting the run, sixteen hex characters
to_text(), from_text()the canonical one-line form, and its strict parser
with_seed, with_hours, with_fault, sampling_every, open_loopderive a variant

The two defaults worth knowing are that a trip ends the run and that the driver does not force IDV(12). Both are Class C quirks of the original that this port fixes by default and can reproduce on request; see the delta register.

Simulation

Simulation(scenario) holds a plant, a controller stack and an integrator state, all of it owned. The original keeps its whole working set in six Fortran COMMON blocks, which allows exactly one simulation per process and no reentrancy at all. This allows as many as there are threads.

run() runs the whole scenario and returns a Run. It runs a copy, so calling it twice gives two equal runs rather than one run and one empty one, and it never raises for a plant that misbehaves: a trip or a failed temperature solve is reported through Run.outcome, because a run that ended early is data.

run() releases the GIL for the entire integration, which is where all the time goes, so an ensemble is a ThreadPoolExecutor and not a ProcessPoolExecutor. Nothing has to be pickled, and the plant does no I/O and touches no Python object while the GIL is down.

Run

MemberWhat it is
to_numpy()the whole run, one (n_samples, 53) float64 array, C-contiguous
hours, stepssimulated time and integrator step per row
measurement(n), manipulated(n)one channel, one-based as XMEAS(n) and XMV(n) are
column(i), columns()one channel by zero-based position, or all 53 keyed by name
labels()ground truth: 'active' and 'since_onset', both (n_samples, 20)
outcome'completed', 'tripped' or 'solve_failed'
tripped_at, tripped_hours, trip_causewhere and why the plant shut down, or None
solve_failed_atthe step a temperature solve failed at, or None; delta D-001
scenariowhat was asked for

Every array is a read-only view over one buffer, which is filled once when the run finishes and moved into NumPy rather than copied. to_numpy() returns the same object every call, and measurement, manipulated, column and columns are strided views into it. Call .copy() for something writable, or numpy.ascontiguousarray if contiguity matters. columns() is a dict in channel order, which makes it a pandas.DataFrame constructor argument as it stands.

labels() is the part with no counterpart in the original. A published Tennessee Eastman dataset is a matrix and a filename, so every detection-delay figure in the literature is computed against whatever onset its author assumed. Here the onset is in the data.

Fault

faults() returns twenty of these. published is the description at teprob.f:172-191 verbatim, five of which say only "Unknown". effect is what the source actually does, which for those five is perfectly explicit: only the physical interpretation was withheld. shape is 'step', 'random' or 'sticking', line names the teprob.f line the fault acts on, and affects_the_plant is False for the three sticking faults, which touch no equation in the model at all. The same table is in The twenty disturbances.

What Scenario(...) cannot say

Three of a scenario's eleven fields are out of the constructors' reach: the event schedule, the continuous extension that allows a fault at a fraction of its full strength, and the choice of integrator. Scenario.from_text is the way to them from Python, and since the text is the scenario's real serialisation, the thing the digest is taken over and the browser app puts in a URL fragment, building through it is building in the form the run will be described by.

text = tep.Scenario.baseline(hours=18).to_text()
scheduled = tep.Scenario.from_text(text.replace("events=", "events=6:start:4,12:stop:4"))

repr() knows about this. A scenario the constructor can express is printed as a constructor call, because that is what is readable at a prompt, and anything else is printed as Scenario.from_text(...). Both round-trip, so eval(repr(s)) == s holds for either shape. The alternative, printing only the constructor's arguments, produced valid Python that evaluated to a different scenario with a different digest and said nothing about it.

notebooks/04-custom-scenarios.ipynb works through all three fields, with a helper that rebuilds the text line from a dictionary rather than patching it, so that a mistyped field name is an error rather than a silently different run.

What is not in the package

tepsim-stats, the crate the repository's own validation ladder runs its statistics on, is development-only. cargo xtask ci asserts that no shipped crate so much as names it, and it is therefore not reachable from Python. That is deliberate: the ladder's statistics have to be free of any dependency the port itself does not have, so they can never be the reason a validation number changed.

For monitoring work the notebooks use notebooks/pcamon.py instead, which is PCA, both monitoring statistics and their control limits in NumPy and the Python standard library alone, with no SciPy and no scikit-learn. Where the two could differ they agree: 02-fault-detection-pca.ipynb checks every printed digit of its detector against the Rust one.

Types

The package ships py.typed, and _tepsim.pyi beside __init__.py describes the compiled half, so mypy and an editor see the real signatures rather than Any.

Worked examples

The four notebooks in notebooks/ are where the package is actually used. They are executed and committed with their outputs and plots, and rendered copies are published beside this book:

Those four links resolve on the published site, where .github/workflows/pages.yml renders the notebooks beside the book with jupyter nbconvert. In a local mdbook build they point at files that are not there; the sources are notebooks/*.ipynb in the repository, and notebooks/README.md says how to run them.

The tutorials in this book are the narrative around them.