18.2 Property-Based Testing and Mutation

1. The Oracle Problem

In 2019, a widely used genomics pipeline reportedly passed all 400 of its example-based tests while silently truncating floating-point precision in a normalization step, producing months of subtly wrong results before anyone noticed. The root cause was not missing tests; it was that every test checked a specific input against a hardcoded expected output, and none asked whether the function obeyed the mathematical properties it was supposed to preserve.

Your test framework just fired ten thousand random inputs at a function and never saw a crash, yet a silent arithmetic error has been corrupting every result, because no assertion knew what the correct output should be. This is the oracle problem: given a random input to a function, how do you determine the right answer? For a sorting function, you can check that the output is sorted and is a permutation of the input without knowing the exact output sequence. For a machine learning model's prediction, there is no simple oracle. The oracle problem explains why most automated test generators produce tests with weak assertions ("does not crash") rather than strong ones ("produces the correct output").

Properties are the solution to the oracle problem. Instead of specifying the exact output for each input, you specify relationships that must hold between the input and output, or between multiple calls to the function. These relationships are called behavioral invariants, and they fall into several families.

What Makes a Good Invariant

A behavioral invariant is a predicate over a function's inputs and outputs that must evaluate to true for every valid input, not just for a handful of chosen examples. Invariants transform testing from a sampling exercise (checking a few points) into a specification exercise (asserting universal laws the implementation must obey). A property-based testing framework like Hypothesis (a Python library for property-based testing that generates random inputs from declarative descriptions called strategies) draws hundreds or thousands of random inputs from a defined domain, runs each through the function, and checks whether the invariant holds. When any input violates the invariant, the framework shrinks the input (progressively simplifies the failing case to find the smallest counterexample) that still triggers the failure. Use behavioral invariants when you cannot write a complete oracle for every input (the common case in scientific and ML code). Fall back to example-based tests only when the expected output is both known and meaningful to assert directly, such as regression snapshots or golden-file comparisons. In short: You do not need to know the right answer to know the answer is wrong; you only need a property it must obey.

"""
The five families of behavioral invariants.
Each family provides a different angle on correctness
without requiring a complete oracle.
"""
from hypothesis import given, settings, assume
from hypothesis import strategies as st
import math


# Family 1: ROUNDTRIP (encode then decode returns the original)
# If f and g are inverses, then g(f(x)) == x for all valid x.

@given(st.floats(min_value=0.01, max_value=1e6, allow_nan=False))
def test_log_exp_roundtrip(x):
    """exp(log(x)) should return x for positive values."""
    assert math.isclose(math.exp(math.log(x)), x, rel_tol=1e-9)


# Family 2: IDEMPOTENCE (applying twice equals applying once)
# f(f(x)) == f(x)

@given(st.lists(st.integers(), min_size=1))
def test_sort_idempotent(xs):
    """Sorting an already-sorted list should produce the same list."""
    once = sorted(xs)
    twice = sorted(once)
    assert once == twice


# Family 3: INVARIANT PRESERVATION (some property holds before and after)
# P(x) implies P(f(x))

@given(st.lists(st.integers()))
def test_sort_preserves_length(xs):
    """Sorting should not change the number of elements."""
    assert len(sorted(xs)) == len(xs)


@given(st.lists(st.integers()))
def test_sort_preserves_elements(xs):
    """Sorting should not add or remove elements."""
    assert sorted(sorted(xs)) == sorted(xs)  # multiset equality


# Family 4: METAMORPHIC (relating outputs of related inputs)
# If input x becomes x', then output f(x) relates to f(x') in a known way.

@given(
    st.lists(st.integers(), min_size=1),
    st.integers()
)
def test_sort_metamorphic_append(xs, new_element):
    """Adding an element and re-sorting should insert it in order."""
    sorted_original = sorted(xs)
    sorted_with_new = sorted(xs + [new_element])
    # The new sorted list should contain all original elements plus the new one
    assert len(sorted_with_new) == len(sorted_original) + 1


# Family 5: ORACLE COMPARISON (comparing against a reference implementation)
# f_optimized(x) == f_reference(x) for all x

def naive_normalize(values):
    """Reference implementation: simple but obviously correct."""
    min_v, max_v = min(values), max(values)
    if max_v == min_v:
        return [0.5] * len(values)
    return [(v - min_v) / (max_v - min_v) for v in values]

def fast_normalize(values):
    """Optimized implementation under test."""
    import numpy as np
    arr = np.array(values, dtype=np.float64)
    mn, mx = arr.min(), arr.max()
    if mx == mn:
        return [0.5] * len(values)
    return ((arr - mn) / (mx - mn)).tolist()

@given(st.lists(
    st.floats(min_value=-1e6, max_value=1e6, allow_nan=False, allow_infinity=False),
    min_size=2
))
def test_normalize_oracle_comparison(values):
    """Optimized normalization matches naive reference."""
    naive = naive_normalize(values)
    fast = fast_normalize(values)
    for n, f in zip(naive, fast):
        assert math.isclose(n, f, rel_tol=1e-9, abs_tol=1e-12)

Once you know which families of properties to test, the next question is how to generate the inputs that will exercise them, especially when those inputs must obey domain-specific constraints like sorted wavelengths or physically valid temperature ranges.

2. Building Hypothesis Strategies for Scientific Data

Hypothesis generates test inputs from strategies: composable descriptions of the input space. Built-in strategies cover primitives (integers, floats, strings, booleans), containers (lists, dicts, tuples), and compositions (mapped, filtered, one-of). For scientific applications, you typically need custom strategies that generate domain-specific data: valid chemical formulas, well-formed molecular graphs, plausible time series, or physically meaningful parameter combinations.

"""
Custom Hypothesis strategies for scientific data types.
"""
from hypothesis import strategies as st
from hypothesis.strategies import composite
from dataclasses import dataclass


@dataclass
class Spectrum:
    """A spectral measurement with wavelengths and intensities."""
    wavelengths: list[float]  # nm, sorted ascending
    intensities: list[float]  # arbitrary units, same length


@composite
def spectra(draw, min_points=2, max_points=200):
    """Generate valid Spectrum objects with sorted wavelengths."""
    n = draw(st.integers(min_value=min_points, max_value=max_points))

    # Wavelengths: sorted, positive, in visible/near-IR range
    wavelengths = sorted(draw(st.lists(
        st.floats(min_value=200.0, max_value=2500.0,
                  allow_nan=False, allow_infinity=False),
        min_size=n, max_size=n
    )))

    # Intensities: non-negative (photon counts)
    intensities = draw(st.lists(
        st.floats(min_value=0.0, max_value=1e6,
                  allow_nan=False, allow_infinity=False),
        min_size=n, max_size=n
    ))

    return Spectrum(wavelengths=wavelengths, intensities=intensities)


@composite
def experiment_configs(draw):
    """Generate plausible experiment configuration dicts."""
    return {
        "temperature": draw(st.floats(
            min_value=4.0, max_value=1500.0  # Kelvin
        )),
        "pressure": draw(st.floats(
            min_value=0.001, max_value=100.0  # atmospheres
        )),
        "concentration": draw(st.floats(
            min_value=1e-9, max_value=1.0  # molar
        )),
        "duration": draw(st.floats(
            min_value=0.1, max_value=86400.0  # seconds (up to 1 day)
        )),
        "replicates": draw(st.integers(
            min_value=1, max_value=96  # standard microplate
        )),
    }


# Using the custom strategies in property tests
@given(spectrum=spectra())
def test_spectrum_wavelengths_sorted(spectrum):
    """Generated spectra should always have sorted wavelengths."""
    for i in range(1, len(spectrum.wavelengths)):
        assert spectrum.wavelengths[i] >= spectrum.wavelengths[i - 1]


@given(spectrum=spectra(), config=experiment_configs())
def test_normalize_preserves_spectrum_length(spectrum, config):
    """Normalization should not change the number of data points."""
    # Assume non-trivial spectrum (not all zeros)
    assume(max(spectrum.intensities) > 0)
    normalized = naive_normalize(spectrum.intensities)
    assert len(normalized) == len(spectrum.intensities)

Property-based testing strengthens the inputs your tests explore, but it cannot tell you whether your assertions are strong enough to catch the faults that matter; for that, you need a way to measure the detection power of the test suite itself.

3. Mutation Testing: Measuring Test Strength

Coverage tells you which code your tests execute. Mutation testing tells you whether your tests would notice if that code were wrong. The process works as follows: make a small change (a mutant) to the source code, then run your test suite. If at least one test fails, the mutant is killed (your tests detected the fault). If all tests pass, the mutant survived (your tests are blind to that class of fault). Figure 18.2 illustrates this feedback loop.

Before we quantify mutation results, one more category of mutant needs definition. An equivalent mutant is one that changes the source code without changing the program's behavior (for example, replacing x + 0 with x). Equivalent mutants are the bane of mutation testing because they inflate the denominator without representing real faults. In practice, most mutation testing tools accept a small rate of equivalent mutants as noise.

The mutation score is the fraction of non-equivalent mutants killed:

3.1 Mutation Operators

A mutation operator defines a class of source-code transformations. Common operators include:

• Arithmetic operator replacement: + becomes -, * becomes /
• Comparison operator replacement: < becomes <=, == becomes !=
• Boolean operator replacement: and becomes or
• Constant replacement: 0 becomes 1, True becomes False
• Statement deletion: removing a line entirely
• Return value mutation: replacing return x with return None

The theoretical foundation is the competent programmer hypothesis: real bugs are small deviations from the correct program. If your tests catch all small deviations (mutants), they will typically also catch most real bugs. The related coupling effect (the empirical observation that tests detecting simple faults also tend to detect more complex faults composed of those simple faults) reinforces this. DeMillo, Lipton, and Sayward (1978) established these hypotheses, and Jia and Harman (2011) validated them empirically, justifying mutation testing as a proxy for real-world fault detection.

"""
Running mutation testing with mutmut on a Python module.
"""
# First, install mutmut: pip install mutmut

# The function under test (spectrum.py):
def find_peaks(
    wavelengths: list[float],
    intensities: list[float],
    threshold: float = 0.1,
    min_distance: int = 5
) -> list[int]:
    """Find peaks in a spectrum above a relative threshold.

    Args:
        wavelengths: Wavelength values (sorted ascending).
        intensities: Intensity values.
        threshold: Minimum relative height (0 to 1) above baseline.
        min_distance: Minimum number of points between peaks.

    Returns:
        Indices of detected peaks.
    """
    if len(intensities) < 3:
        return []

    baseline = min(intensities)
    max_intensity = max(intensities)
    abs_threshold = baseline + threshold * (max_intensity - baseline)

    peaks = []
    for i in range(1, len(intensities) - 1):
        # A peak is a local maximum above the threshold
        if (intensities[i] > intensities[i - 1] and
                intensities[i] > intensities[i + 1] and
                intensities[i] >= abs_threshold):
            # Enforce minimum distance from the last peak
            if not peaks or (i - peaks[-1]) >= min_distance:
                peaks.append(i)

    return peaks

With the function under test defined, we can use mutmut (a Python mutation testing tool that automatically applies mutation operators to source code and re-runs the test suite against each mutant) to measure how many of those decision points our tests actually guard.

# Run mutmut on the spectrum module
mutmut run --paths-to-mutate=spectrum.py --tests-dir=tests/

# View results summary
mutmut results

# Example output:
# Legend for output:
# -Ass = survived (test suite did not catch the mutant)
# - Killed = killed (test suite caught the mutant)
# - Timeout = timeout (mutant caused an infinite loop)
#
# Survived:   3
# Killed:    24
# Timeout:    1
# Total:     28
#
# Mutation score: 24/27 = 88.9% (excluding timeouts)

# Inspect surviving mutants
mutmut show 5
# --- spectrum.py
# +++ spectrum.py (mutant 5)
# @@ -15,7 +15,7 @@
# -        if (intensities[i] > intensities[i - 1] and
# +        if (intensities[i] >= intensities[i - 1] and
#
# This mutant changed > to >= in the local maximum check.
# It survives because no test has a plateau (equal adjacent values)
# where the distinction between > and >= matters.

Each surviving mutant points to a gap in the test suite. The mutant above (changing > to >= in the local maximum check) survives because no test includes a spectrum with a plateau. This is precisely the kind of edge case that property-based testing excels at finding: flat regions, equal adjacent values, monotonic sequences.

3.2 Killing Surviving Mutants with Property Tests

Mutation testing and property-based testing reinforce each other: mutation testing identifies where the suite is weak, and property-based testing generates inputs that exercise those weak spots. The workflow is: run mutmut, examine surviving mutants, write property tests targeting the uncovered behavior, and re-run to confirm the kills.

"""
Killing surviving mutants with targeted property tests.
"""
from hypothesis import given, assume
from hypothesis import strategies as st


# Mutant 5 survived: > changed to >= in local max check.
# This means no test has equal adjacent values where a "peak"
# would be incorrectly detected. Write a property test for plateaus.

@given(st.lists(
    st.floats(min_value=0.0, max_value=100.0,
              allow_nan=False, allow_infinity=False),
    min_size=5, max_size=50
))
def test_plateau_not_detected_as_peak(intensities):
    """A flat region (constant value) should never contain a peak.

    A peak requires being strictly greater than both neighbors.
    If we make a subarray constant, no index in that subarray
    should appear in the peak list.
    """
    # Create a spectrum with a flat plateau in the middle
    n = len(intensities)
    mid = n // 2
    plateau_value = intensities[mid]
    plateau_start = max(1, mid - 2)
    plateau_end = min(n - 1, mid + 3)

    # Set plateau region to a constant value
    for i in range(plateau_start, plateau_end):
        intensities[i] = plateau_value

    wavelengths = list(range(n))  # dummy wavelengths
    peaks = find_peaks(wavelengths, intensities, threshold=0.0)

    # No peak should be in the plateau interior
    for p in peaks:
        if plateau_start < p < plateau_end - 1:
            # This is interior: neighbors on both sides have the same value
            assert intensities[p] > intensities[p - 1] or \
                   intensities[p] > intensities[p + 1], \
                f"Peak at index {p} is on a plateau (value={intensities[p]})"


# Mutant 12 survived: min_distance changed from >= to >
# This means no test has peaks exactly min_distance apart.

@given(
    st.integers(min_value=3, max_value=10),
    st.floats(min_value=1.0, max_value=100.0,
              allow_nan=False, allow_infinity=False)
)
def test_peaks_at_exact_min_distance(min_dist, peak_height):
    """Two peaks exactly min_distance apart should both be detected."""
    # Construct a spectrum with two peaks at exactly min_distance spacing
    n = min_dist * 3 + 2
    intensities = [0.0] * n
    peak1_idx = min_dist
    peak2_idx = peak1_idx + min_dist  # exactly min_distance apart

    intensities[peak1_idx] = peak_height
    intensities[peak2_idx] = peak_height

    wavelengths = list(range(n))
    peaks = find_peaks(
        wavelengths, intensities,
        threshold=0.0, min_distance=min_dist
    )

    assert peak1_idx in peaks, f"First peak at {peak1_idx} not detected"
    assert peak2_idx in peaks, f"Second peak at {peak2_idx} not detected"

Key Insight: Mutation Score Measures Assertion Strength, Not Coverage

A test suite can have 100% branch coverage and a 60% mutation score. That means 40% of the injected faults go undetected: the tests reach the right code but do not check the right properties. Mutation score is a strictly stronger quality metric than coverage because it requires not just execution but detection. For scientific software where an undetected arithmetic error could invalidate published results, mutation score is in practice the metric that matters. A common guideline is to aim for at least 80% mutation score on critical modules; scores below 70% typically indicate that the test suite is fundamentally incomplete.

4. Combining Property Tests and Mutation Testing in Continuous Integration (CI)

In a CI pipeline, property tests and mutation tests serve complementary roles. Property tests run on every commit (they are fast, typically under a minute for 100 properties with 100 examples each). Mutation tests run nightly or on pull requests (they are slow, requiring one full test suite run per mutant). The combination provides both rapid feedback (property tests catch regressions immediately) and deep analysis (mutation tests verify that the test suite itself remains strong).

"""
CI configuration for combined property and mutation testing.
This pytest conftest.py configures Hypothesis for CI environments.
"""
from hypothesis import settings, Phase, HealthCheck

# CI profile: more examples, deterministic, no shrinking timeout
settings.register_profile(
    "ci",
    max_examples=500,          # more thorough than dev (default 100)
    derandomize=True,          # reproducible failures
    suppress_health_check=[
        HealthCheck.too_slow,  # allow slow strategies in CI
    ],
    phases=[
        Phase.explicit,        # run explicit examples first
        Phase.reuse,           # replay previously failing examples
        Phase.generate,        # generate new examples
        Phase.shrink,          # shrink failing examples to minimal form
    ],
    deadline=None,             # no per-example time limit
)

# Dev profile: fast feedback
settings.register_profile(
    "dev",
    max_examples=20,           # fast iteration
    derandomize=False,         # explore new inputs each run
    deadline=2000,             # 2-second timeout per example
)

# Load profile from environment variable:
#   HYPOTHESIS_PROFILE=ci pytest tests/
settings.load_profile("dev")  # default to dev

# .github/workflows/test.yml (excerpt)
# Property tests on every push; mutation tests on PRs only.

jobs:
  property-tests:
    runs-on: ubuntu-latest
    steps:
      - uses: actions/checkout@v4
      - uses: actions/setup-python@v5
        with:
          python-version: "3.12"
      - run: pip install -r requirements-test.txt
      - run: HYPOTHESIS_PROFILE=ci pytest tests/ -x --tb=short

  mutation-tests:
    if: github.event_name == 'pull_request'
    runs-on: ubuntu-latest
    steps:
      - uses: actions/checkout@v4
      - uses: actions/setup-python@v5
        with:
          python-version: "3.12"
      - run: pip install -r requirements-test.txt mutmut
      - run: |
          mutmut run \
            --paths-to-mutate=src/ \
            --tests-dir=tests/ \
            --runner="pytest tests/ -x --tb=no -q"
      - run: |
          # Fail the PR if mutation score drops below 80%
          python -c "
          import json
          with open('.mutmut-cache/results.json') as f:
              r = json.load(f)
          killed = r.get('killed', 0)
          total = r.get('total', 1) - r.get('equivalent', 0)
          score = killed / total * 100 if total > 0 else 100
          print(f'Mutation score: {score:.1f}%')
          assert score >= 80, f'Mutation score {score:.1f}% below 80% threshold'
          "

5. The Metamorphic Testing Pattern for Machine Learning (ML) Systems

ML models present a particularly severe oracle problem: for a given input, no formula can tell you the "correct" prediction. Metamorphic testing (introduced as Family 4 in the invariant taxonomy above) sidesteps this by testing relationships between predictions rather than individual predictions. If rotating an image by 5 degrees changes a classifier's prediction from "cat" to "dog," that is almost certainly a bug, even though we do not know the correct label. This section develops the metamorphic pattern in depth, showing how to choose relations, compose them with Hypothesis, and interpret failures, because metamorphic testing is the primary technique for validating the ML-intensive systems we build in Part III and the evaluation frameworks in Chapter 56.

"""
Metamorphic testing for a scientific classifier.
Tests relationships between predictions, not individual predictions.
"""
from hypothesis import given, assume
from hypothesis import strategies as st
import numpy as np


def classify_material(spectrum: np.ndarray) -> str:
    """Classify a material from its infrared (IR) absorption spectrum.
    (Placeholder for an actual ML model.)
    """
    # In production, this calls model.predict()
    peak_position = np.argmax(spectrum)
    if peak_position < len(spectrum) * 0.3:
        return "polymer"
    elif peak_position < len(spectrum) * 0.7:
        return "ceramic"
    else:
        return "metal"


# Metamorphic relation 1: Scaling invariance
# Multiplying all intensities by a positive constant should not
# change the classification (the material is the same, just measured
# with a different detector gain).

@given(
    st.lists(
        st.floats(min_value=0.1, max_value=1000.0,
                  allow_nan=False, allow_infinity=False),
        min_size=20, max_size=200
    ),
    st.floats(min_value=0.01, max_value=100.0,
              allow_nan=False, allow_infinity=False)
)
def test_classification_scale_invariant(intensities, scale_factor):
    """Scaling intensities should not change the classification."""
    original = np.array(intensities)
    scaled = original * scale_factor

    assert classify_material(original) == classify_material(scaled), \
        f"Classification changed when scaling by {scale_factor}"


# Metamorphic relation 2: Noise robustness
# Adding small Gaussian noise should not change the classification
# (real measurements always have some noise).

@given(
    st.lists(
        st.floats(min_value=1.0, max_value=1000.0,
                  allow_nan=False, allow_infinity=False),
        min_size=20, max_size=200
    ),
    st.integers(min_value=0, max_value=2**31 - 1)  # random seed
)
def test_classification_noise_robust(intensities, seed):
    """Small noise should not change the classification."""
    rng = np.random.default_rng(seed)
    original = np.array(intensities)
    snr = 100.0  # signal-to-noise ratio
    noise = rng.normal(0, original.max() / snr, size=len(original))
    noisy = original + noise

    assert classify_material(original) == classify_material(noisy), \
        f"Classification changed with SNR={snr}"

Each technique covered so far (property tests, mutation analysis, metamorphic relations) captures a different dimension of test quality, raising the question of how to unify these signals into a single actionable metric.

6. From Test Metrics to Test Confidence

Three complementary metrics capture test quality: coverage (what code executes), mutation score (what faults are detected), and property count (how many invariants are verified). None alone is sufficient. A useful composite metric is the mutation-adjusted confidence score:

where \(C_{\text{branch}}\) is branch coverage (0 to 1), \(M\) is mutation score (0 to 1), and \(P / P_{\text{target}}\) is the ratio of property tests to a target count (capped at 1). A module with 95% branch coverage, 90% mutation score, and 8 out of 10 target properties gets a confidence of \(0.95 \times 0.90 \times 0.80 = 0.684\). This composite score makes the tradeoffs visible: high coverage with low mutation score (weak assertions) pulls the score down just as much as low coverage with strong assertions.

Section 18.3 builds this metric into an automated invariant discovery system that targets all three components simultaneously.