Chapter 4 — Baseline replication and pilot

Lesson map

Lesson 04 published. Map, example, EX021.

flowchart TD
  N0["Learning outcomes and connection"]
  N1["4.1 Reproduction, replication, robustness"]
  N2["4.2 Reproducing a model calculation: the inoperabili"]
  N3["4.3 The GridWater-Sim discrepancy synthetic"]
  N4["4.4 Verifying a simulator"]
  N5["4.5 Replicating a detection baseline: evaluation art"]
  N6["4.6 Eliciting parameters"]
  N7["4.7 The pilot"]
  N0 --> N1 --> N2 --> N3 --> N4 --> N5 --> N6 --> N7
  N7 --> EX["Worked example"]
  EX --> Q["Exercise EX021"]

Logical drawing

flowchart LR
  N0["Learning outcomes and connection"]
  N1["4.1 Reproduction, replication, robustness"]
  N2["4.2 Reproducing a model calculation: the inoperabili"]
  N3["4.3 The GridWater-Sim discrepancy synthetic"]
  N4["4.4 Verifying a simulator"]
  N5["4.5 Replicating a detection baseline: evaluation art"]
  N6["4.6 Eliciting parameters"]
  N7["4.7 The pilot"]
  N0 --> N1 --> N2 --> N3 --> N4 --> N5 --> N6 --> N7 --> Q["EX021"]

[Learning outcomes and connection]
    |
    v
[4.1 Reproduction, replication, robustness]
    |
    v
[4.2 Reproducing a model calculation: the inoperabili]
    |
    v
[4.3 The GridWater-Sim discrepancy synthetic]
    |
    v
[4.4 Verifying a simulator]
    |
    v
[4.5 Replicating a detection baseline: evaluation art]
    |
    v
[4.6 Eliciting parameters]
    |
    v
[4.7 The pilot]
    |
    v
[Worked example] --> [Exercise EX021]

Learning outcomes and connection

By the end of this chapter the learner will be able to:

1. Distinguish reproduction, replication and robustness checks, and choose which one a baseline needs.
2. Reproduce a published model calculation, predeclare a tolerance, diagnose a discrepancy and write a reproduction log.
3. Verify a simulator with invariants, unit tests and an independent implementation before trusting any comparison it produces.
4. Recognize evaluation-protocol artifacts — point adjustment, temporal leakage, testbed-specific features — that inflate detection baselines.
5. Elicit and document model parameters from experts with explicit uncertainty.
6. Design a bounded pilot whose purpose is feasibility, with go/no-go criteria, and explain why a pilot must not be used to estimate the effect size.

Chapter 3 fixed the hypotheses, baselines and failure conditions. Before the main study can compare DEP with SILO, CRITFLAG and ORACLE, every number those comparisons depend on must be shown to mean what it claims.

4.1 Reproduction, replication, robustness

The National Academies (2019) distinguish reproducibility — obtaining consistent results using the same input data, computational steps, methods and conditions of analysis — from replicability — obtaining consistent results across studies aimed at the same question, each with its own data. A third activity, often called a robustness check, varies analytic choices on the same data. Infrastructure theses need all three:

Lesson table; scroll horizontally to read all columns →

4.2 Reproducing a model calculation: the inoperability input–output model

4.2.1 The model

Haimes and Jiang's inoperability input–output model (2001) adapts Leontief's input–output economics to infrastructure. Each sector i has an inoperability q_i between 0 (fully operable) and 1 (fully inoperable). The equilibrium after an initial perturbation c* is q = A*q + c*, hence q = (I − A*)⁻¹ c*. The element A*_ij is the inoperability contributed to sector i per unit inoperability of sector j. The convention matters: the row is the receiving sector, the column the supplying sector. Papers and software differ in notation, and some derive A* from economic input–output tables through normalization steps that are easy to get subtly wrong.

4.2.2 Reproduction lab

A synthetic "published" paper reports that, after a perturbation making power 50% inoperable, water reaches an equilibrium inoperability of 0.236 in a three-sector model. Predeclare a tolerance of ± 0.005 (rounding in the publication), then reproduce:

# Lab 4a - reproducing an inoperability input-output (IIM) result (offline, standard library only)
# q = (I - A*)^-1 c*. Convention: A*[i][j] = inoperability contributed to sector i per unit inoperability of sector j.
def solve(a, c):
    n = len(c)
    m = [[(1.0 if i == j else 0.0) - a[i][j] for j in range(n)] + [c[i]] for i in range(n)]
    for col in range(n):
        piv = max(range(col, n), key=lambda r: abs(m[r][col]))
        m[col], m[piv] = m[piv], m[col]
        for r in range(n):
            if r != col:
                f = m[r][col] / m[col][col]
                m[r] = [x - f * y for x, y in zip(m[r], m[col])]
    return [m[i][n] / m[i][i] for i in range(n)]
sectors = ["power", "water", "telecom"]
A = [[0.00, 0.05, 0.10],   # power depends a little on water (cooling) and telecom (SCADA comms)
     [0.40, 0.00, 0.15],   # water depends strongly on power
     [0.30, 0.02, 0.00]]   # telecom depends on power
c = [0.50, 0.00, 0.00]    # initial perturbation: power 50% inoperable
q = solve(A, c)
print("convention A* (row = receiving sector):", [round(x, 3) for x in q])
At = [list(r) for r in zip(*A)]
qt = solve(At, c)
print("transposed by mistake:               ", [round(x, 3) for x in qt])

Expected output:

convention A* (row = receiving sector): [0.528, 0.236, 0.163]
transposed by mistake:                [0.528, 0.028, 0.057]

With the correct convention the reproduction matches (water 0.236). With the matrix transposed — an error that produces no warning, no crash and a plausible-looking vector — water's inoperability falls by almost an order of magnitude, and the ranking of sectors changes: telecom now appears more affected than water. Power's own value happens to agree in both cases, which is exactly what makes the error dangerous: a replicator who checks only the perturbed sector would declare success.

4.2.3 The reproduction log

A reproduction log records: the target value and its source (page, table); the tolerance, fixed before running; software environment (language version, operating system); inputs and their provenance; the exact output; discrepancies and their diagnosis; and the decision (reproduced, reproduced after correction, not reproduced). When the original authors' convention is ambiguous, write to them; record the correspondence. A non-reproduction is a result, not an embarrassment, and belongs in the thesis.

4.2.4 What a reproduced IIM can and cannot support

Reproducing the calculation shows only that the arithmetic is understood. The IIM remains linear and equilibrium-based: it does not model buffers, timing or restoration order, and its coefficients are usually derived from economic flows rather than engineering dependencies. In the Cascade Study it serves as a sector-level cross-check — do the sectors the IIM ranks as most affected match the services the buffer model shows failing? — not as a competitor to the restoration model.

4.3 The GridWater-Sim discrepancy (synthetic)

The synthetic benchmark GridWater-Sim reports, for a Riverbend-like configuration, a mean weighted critical-service loss under customer-first restoration that the candidate's re-implementation could not match within the predeclared tolerance. Diagnosis followed the order that most often finds the cause:

1. Inputs. Were the same parameter distributions used? The benchmark documentation gave triangular distributions as (low, mode, high); Python's random.triangular takes (low, high, mode). Swapping mode and high silently shifts the mean.
2. Semantics. Does "restored" mean re-energized or service resumed? GridWater-Sim applied a one-hour restart lag only to water; the re-implementation applied it to every service.
3. Numerics. The benchmark advanced time in fixed steps; the re-implementation computed exact outage intervals. Fixed steps introduce discretization error that depends on the step and on where buffer expiries fall within a step.
4. Randomness. Were seeds and draw order identical? Draw order matters: drawing durations before buffers in one implementation and after in another produces different scenario sets from the same seed.

After aligning (1) and (2), the remaining difference was within tolerance and shrank with smaller time steps, confirming (3). The thesis records the benchmark as reproduced after correcting a parameter-order ambiguity, reports the original and corrected values, and adopts the exact-interval implementation with the time-step version kept as an independent check.

4.4 Verifying a simulator

Verification asks whether the model is implemented correctly; validation asks whether it represents the real system adequately for its purpose (Sargent, 2013). Verification comes first, because validating a buggy implementation is meaningless.

Invariants that any correct restoration simulator must satisfy:

• If every buffer exceeds the horizon, weighted critical-service loss is zero under every ordering.
• Two policies that produce the same order produce identical results on identical draws.
• The oracle's loss is less than or equal to every policy's loss on every draw.
• Increasing a buffer, holding everything else fixed, never increases loss.
• Customer-hours under Smith's rule are less than or equal to those under any other ordering when restoration durations are equal.
• Loss and customer-hours are non-negative and bounded by weight × horizon and customers × horizon.

Unit tests check the smallest components: interval arithmetic for a single service; the effect of a dependency with an on-site buffer; crew assignment with two crews. Cross-implementation tests compare two independently written implementations on fixed draws. Extreme-condition tests use zero buffers, zero durations and one feeder. Keep the tests in the archive (Chapter 7); a reviewer who can run them trusts the comparisons more than one who must take them on faith.

4.5 Replicating a detection baseline: evaluation artifacts

Even though the Cascade Study does not propose a detector, its time-to-impact analysis depends on when a disruption is recognized, and the candidate must be able to judge detection baselines reported for infrastructure. The common inflation mechanisms are evaluation choices, not model properties.

4.5.1 Point adjustment

Kim et al. (2022) showed that point adjustment (PA) — counting a whole anomalous segment as detected when any point in it is flagged — can make random scores look excellent. The following lab reproduces the effect on synthetic labels:

# Lab 4b - how point adjustment (PA) inflates detection scores (offline, standard library only)
import random
random.seed(42)
N = 2000
labels = [0] * N
for start, length in [(300, 60), (900, 120), (1500, 80)]:
    for t in range(start, start + length):
        labels[t] = 1
scores = [random.random() for _ in range(N)]          # a detector that knows nothing
thr = sorted(scores)[int(0.97 * N)]                    # alarm on the top 3% of scores
pred = [1 if s > thr else 0 for s in scores]
def f1(p, y):
    tp = sum(1 for a, b in zip(p, y) if a and b)
    fp = sum(1 for a, b in zip(p, y) if a and not b)
    fn = sum(1 for a, b in zip(p, y) if b and not a)
    prec = tp / (tp + fp) if tp + fp else 0.0
    rec = tp / (tp + fn) if tp + fn else 0.0
    return round(2 * prec * rec / (prec + rec), 3) if prec + rec else 0.0
def segments(y):
    segs, t = [], 0
    while t < len(y):
        if y[t]:
            s = t
            while t < len(y) and y[t]:
                t += 1
            segs.append((s, t))
        else:
            t += 1
    return segs
def adjust(p, y, k=0.0):
    q = list(p)
    for s, e in segments(y):
        hit = sum(p[s:e]) / (e - s)
        if hit > k:
            for t in range(s, e):
                q[t] = 1
    return q
print("point-wise F1           :", f1(pred, labels))
print("F1 after PA             :", f1(adjust(pred, labels), labels))
print("F1 after PA%K (K = 0.20):", f1(adjust(pred, labels, 0.20), labels))
print("segments hit by >=1 alarm:", sum(1 for s, e in segments(labels) if any(pred[s:e])), "of", len(segments(labels)))

Expected output:

point-wise F1           : 0.05
F1 after PA             : 0.911
F1 after PA%K (K = 0.20): 0.05
segments hit by >=1 alarm: 3 of 3

A detector with no information scores 0.911 after PA because each long segment contains, by chance, at least one of the 3% random alarms. PA%K with K = 20% removes the inflation here. For CI research the lesson goes further: neither metric says when the first true alarm occurred relative to the segment's start, which is what matters for buffers. A consequence-relevant metric is detection latency relative to time-to-impact — the fraction of incidents detected before the earliest dependent service's buffer is exhausted.

4.5.2 Other artifacts

• Temporal leakage. Random train/test splits of time series let a model learn from the future. Use time-ordered splits that mimic deployment.
• Testbed-specific features. A detector may learn signatures of the testbed's attack-injection tooling or its specific sensors rather than of attacks.
• Attack-set coverage. A dataset's attacks were chosen by its builders; performance on them says little about other attack classes (Arp et al., 2022).
• Base rates. A false-positive rate that looks small on a balanced dataset can swamp operators at realistic base rates (Sommer and Paxson, 2010).

4.6 Eliciting parameters

The model's buffers, restoration durations and weights come from people. Structured expert elicitation — preparing experts, asking for quantiles rather than single values, feeding back implied distributions and recording rationales — reduces anchoring and overconfidence (O'Hagan et al., 2006). For the Cascade Study (synthetic):

Lesson table; scroll horizontally to read all columns →

Triangular distributions are a pragmatic choice for three-point elicitations; they are not physically derived. Chapter 6 tests whether conclusions depend on the distribution family.

4.7 The pilot

4.7.1 What a pilot is for

A pilot tests feasibility: whether procedures, materials, timing, data flows and instruments work. It is not a small main study. Effect sizes estimated from pilots are too imprecise to plan sample sizes or to support claims (Leon, Davis and Kraemer, 2011). The Cascade Study pilot had two strands:

1. Simulation pilot. Run all four policies on 50 draws to check that the code, logging and invariants work, that runtime allows the planned number of draws, and that censoring at 72 hours is rare. The pilot results (differences between policies) are logged as exploratory and not interpreted.
2. Tabletop pilot. Three operators (synthetic participants P01–P03, not eligible for the main study) worked through a 40-minute scenario with and without a time-to-impact display. Questions: Do they understand the display? Does the scenario reveal the outcome too early? How long does each decision take? Does the consent and debrief process work?

4.7.2 Go/no-go criteria (predeclared)

Lesson table; scroll horizontally to read all columns →

Every revision is logged in the deviation log with the date and confirmation that no main-study data existed.

4.8 Worked example: the pilot that changed the scenario

In the tabletop pilot, P02 asked why the hospital's generator buffer was shown as 26 hours when "everyone knows the water goes first". The candidate had displayed only direct power buffers. The display was revised to show effective time-to-impact (the minimum over all paths, as in Lab 3) alongside direct buffers, with the path shown. This changed the instrument — a legitimate pilot outcome — and raised a methodological point for the main study: the comparison between "with display" and "without display" now tested the effective-time-to-impact display, not merely any criticality information. The hypotheses for the tabletop were refined accordingly, before the main study, and the change was preregistered as an amendment.

4.9 Common errors and limits

• Declaring a reproduction successful after checking a single output.
• Setting tolerances after seeing the discrepancy.
• Treating the benchmark's code as the specification when it disagrees with its paper, or vice versa, without recording which was followed.
• Validating a simulator before verifying it.
• Reporting point-adjusted F1 for infrastructure detection without latency relative to consequence.
• Asking experts for single numbers and treating them as certain.

• Using pilot effect sizes to plan or to claim.

• Checking only the perturbed sector in a reproduction, as with power in the transposed inoperability matrix, which agrees under both conventions.

• Assuming a library's parameter order matches the paper's notation, as with the triangular distribution in the GridWater-Sim case.
• Interpreting differences between policies in a pilot, when the pilot's purpose is feasibility and its results are logged as exploratory.
• Changing a tabletop instrument after the pilot without recording the change in the deviation log.

Exercises — chapter 4

EX021 — Reproduction log. Objective: document a reproduction attempt. Inputs: Lab 4a. Constraints: run offline; tolerance ± 0.005 fixed before running. Steps: (1) record the target (water 0.236); (2) run the transposed matrix first and record the discrepancy; (3) diagnose using the matrix convention; (4) run the correct convention; (5) write the log entry with environment details. Expected evidence: reproduction log. Rubric: 3 = exact outputs, convention diagnosis, environment recorded, notes that power agrees in both runs; 2 = one element missing; 1 = no diagnosis; 0 = not run.

EX022 — IIM perturbation analysis. Objective: use a reproduced model for a bounded claim. Inputs: Lab 4a with the correct convention. Constraints: offline. Steps: (1) compute q for perturbations power 0.5 and telecom 0.5 separately; (2) compute q for both together; (3) check whether the joint result equals the sum of the separate results and explain why; (4) state one claim the model supports and one it cannot. Expected evidence: three vectors and a paragraph. Rubric: 3 = correct vectors, linearity (superposition) explained, claim limited to equilibrium inoperability without timing; 2 = vectors correct, explanation thin; 1 = errors; 0 = none.

EX023 — Parameter-order trap. Objective: detect a silent input error. Inputs: a distribution documented as (low = 4, mode = 8, high = 12). Constraints: use Python's random.triangular(low, high, mode). Steps: (1) compute the theoretical mean, (low + mode + high) / 3, for the correct call triangular(4, 12, 8); (2) estimate the mean of the mistaken call triangular(4, 8, 12) from 100,000 seeded draws and record the minimum and maximum; explain what the mistaken call actually samples; state which test would have caught it. Expected evidence: two means and explanation. Rubric: 3 = correct mean 8.0; the mistaken call treats 8 as high and 12 as mode, so CPython's formula returns 4 + 4·√(2u) for uniform u — values in [4, ≈ 9.66] with mean ≈ 7.77, not a triangular distribution at all; a mean check alone is weak, while checking that the sample maximum approaches 12 and plotting a histogram catches it; 2 = error identified, test weak; 1 = error missed; 0 = none.

EX024 — Invariant tests. Objective: write verification tests. Inputs: the invariants in §4.4 and the simulator you will meet in Chapter 5 (Lab 5), or your own. Constraints: at least four invariants as executable assertions. Steps: (1) choose at least four invariants from §4.4; (2) write each as an executable assertion; (3) run the assertions on fixed draws; (4) deliberately introduce a bug, such as removing the restart lag for one service, and show which test fails. Expected evidence: test code and output before and after the bug. Rubric: 3 = four or more meaningful invariants, the bug caught, explanation of why; 2 = tests run but the bug is not caught; 1 = trivial tests; 0 = none.

EX025 — Latency-relative detection metric. Objective: define a consequence-relevant detection metric. Inputs: Lab 4b and Lab 3's time-to-impact values (the earliest critical impact at 2 h for water SCADA, 8 h for telecom). Constraints: define the metric formally. Steps: (1) define "detected before impact" for a segment with known start and a known time-to-impact; (2) compute it for the random detector, treating each time step as 6 minutes and taking time-to-impact as 2 h for every segment; (3) compare with the F1 values. Expected evidence: definition, computed fraction, paragraph. Rubric: 3 = sound definition and correct computation from the first alarm index in each segment, with the contrast to PA explained; 2 = definition fine, computation error; 1 = vague; 0 = none.

EX026 — Pilot plan. Objective: design a feasibility pilot. Inputs: your own study from earlier exercises. Constraints: at least six go/no-go criteria with thresholds; explicit statement of what will not be estimated. Steps: (1) list the pilot's feasibility questions; (2) set predeclared thresholds; (3) define what happens on "no-go" for each; (4) define who may participate and why pilot participants are excluded from the main study. Expected evidence: pilot plan. Rubric: 3 = criteria measurable, feasibility-only, exclusion justified, revision process logged; 2 = one criterion vague; 1 = pilot used to estimate effect; 0 = none.

Chapter synthesis

Before comparing methods, make sure the numbers mean what they claim. Reproduce baselines with a tolerance fixed in advance, check more than one output, and record conventions — a transposed matrix produced a plausible, wrong answer that agreed on the one value a hurried replicator would check. Diagnose discrepancies through inputs, semantics, numerics and randomness. Verify simulators with invariants, unit tests and an independent implementation before validating them. Treat detection baselines with suspicion when they rely on point adjustment, random splits or a single testbed, and prefer latency relative to time-to-impact for infrastructure. Elicit parameters with uncertainty and provenance. Use the pilot to fix procedures and instruments, never to estimate effects.