From 75a949ac1f6559c36a27a3c1a044678ca44a7443 Mon Sep 17 00:00:00 2001
From: Julian Christian Sanders <209085623+JROChub@users.noreply.github.com>
Date: Wed, 9 Sep 2026 16:41:00 -0700
Subject: [PATCH] Add first-divergence diagnostics for particle-track
 comparisons

Compare recorded histories through the existing Tracks reader, with exact
floating-point bits and finite-value ULP distances. Align numeric history
identifiers and report particle/state structure separately from values.

Validate HDF5 schemas and offsets before reading particle segments, retain
track 3.0 and 3.1 compatibility, and distinguish invalid input from mismatch.
Add focused tests and developer usage documentation without changing
transport or its random-number stream.
---
 docs/source/devguide/tests.rst          |  80 ++++
 tests/unit_tests/test_compare_tracks.py | 525 ++++++++++++++++++++++++
 tools/dev/compare_tracks.py             | 231 +++++++++++
 3 files changed, 836 insertions(+)
 create mode 100644 tests/unit_tests/test_compare_tracks.py
 create mode 100644 tools/dev/compare_tracks.py

diff --git a/docs/source/devguide/tests.rst b/docs/source/devguide/tests.rst
index 8627cbde232..0f0208ac9f8 100644
--- a/docs/source/devguide/tests.rst
+++ b/docs/source/devguide/tests.rst
@@ -89,6 +89,86 @@ that, consider the following:
   limit the number of threads that OpenBLAS uses internally; this can be done by
   setting the :envvar:`OPENBLAS_NUM_THREADS` environment variable to 1.
 
+Comparing Recorded Particle Histories
+------------------------------------
+
+``tools/dev/compare_tracks.py`` compares two completed, self-contained version-3.0
+or 3.1 track files using :class:`openmc.Tracks`. It reports the first difference
+in the recorded states without rerunning transport or changing its random-number
+stream. From a checkout with the OpenMC Python package installed, run::
+
+    python tools/dev/compare_tracks.py reference/tracks.h5 candidate/tracks.h5
+
+The JSON report has a ``status`` of ``match``, ``mismatch``, or ``invalid_input``.
+The corresponding exit codes are 0, 1, and 2. Save the report using shell
+redirection if needed. Both input files are read-only. Invalid file types,
+unsupported versions, malformed state schemas, and inconsistent particle
+offsets are reported as invalid input, rather than treated as matching data.
+Two files with no recorded histories are also invalid for this comparison.
+Version-3.0 particle indices are normalized to the corresponding PDG numbers
+used in version 3.1, so equivalent particle types compare equal across formats.
+
+Histories are aligned by their numeric ``(batch, generation, particle number)``
+identifier, independent of HDF5 dataset order. Within a history, particle tracks
+and their states are compared in recorded order. State fields are compared in
+the order ``r.x``, ``r.y``, ``r.z``, ``u.x``, ``u.y``, ``u.z``, ``E``, ``time``,
+``wgt``, ``cell_id``, ``cell_instance``, and ``material_id``. The earliest
+differing state takes precedence over field order. Common prefixes are checked
+before reporting an extra particle track or state. Missing histories, differing
+particle types, and differing record counts have distinct report kinds.
+
+A numerical mismatch includes the history identifier, zero-based
+``particle_index`` and ``state_index``, field name, values from both runs, and
+the binary64 bit patterns for floating-point fields. Finite floating-point
+values also include ``ulp_distance``: the number of representable steps between
+them, counting the two signed zeros as the same numeric point. For example,
+``1.0`` and its next larger representable value differ by one step. Signed zeros
+nevertheless compare unequal because their bits differ. NaNs and infinities
+are compared by bits, including NaN payloads, and have no finite ULP distance
+(``null`` in JSON). Floating-point values are represented as strings so that
+nonfinite values remain valid JSON. Integer fields compare exactly.
+
+The comparison ignores HDF5 padding, byte order, and incidental file attributes
+such as timestamps. It does not compare complete HDF5 file bytes. The reader
+loads both files into memory; select a small set of histories for a focused
+diagnostic run.
+
+For a cross-platform investigation, start with serial, single-threaded
+fixed-source runs and explicit track identifiers. For example, configure the
+same model before exporting it on each system:
+
+.. code-block:: python
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 2
+    model.settings.particles = 100
+    model.settings.seed = 1
+    model.settings.event_based = False
+    model.settings.shared_secondary_bank = False
+    model.settings.track = [(1, 1, 1), (1, 1, 2), (2, 1, 1)]
+
+Use the same source distribution and nuclear data, build with
+``OPENMC_ENABLE_STRICT_FP=ON``, and set ``OMP_NUM_THREADS=1``. Explicit
+:attr:`openmc.Settings.track` identifiers avoid relying on which histories a
+parallel run encounters before reaching ``max_tracks``. Repeat each native
+configuration before comparing platforms. Retain the source revision, compiler
+commands, build configuration, executable and input digests, dependency
+versions, RNG settings, and floating-point environment with each result.
+Strict compiler settings do not themselves make different system math libraries
+produce identical values.
+
+With the local secondary bank used above, ``particle_index=0`` denotes the
+primary particle and subsequent indices denote the recorded secondary
+execution sequence. These indices are not persistent genealogy IDs. Shared
+secondary-bank transport uses separate particle identifiers and should be
+investigated separately. The tool reports the first *recorded-state*
+difference, not necessarily the first differing numerical operation: track
+files do not contain operation operands, reaction identifiers, or RNG state.
+Use the reported history and its preceding matching states to narrow further
+debugger or operation-level tracing without introducing additional RNG draws.
+A match establishes equality of the selected recorded histories, not of every
+unrecorded operation or tally in the calculation.
+
 Debugging Tests in CI
 ---------------------
 
diff --git a/tests/unit_tests/test_compare_tracks.py b/tests/unit_tests/test_compare_tracks.py
new file mode 100644
index 00000000000..3f55752449d
--- /dev/null
+++ b/tests/unit_tests/test_compare_tracks.py
@@ -0,0 +1,525 @@
+"""Exact track diagnostics using self-contained, synthetic HDF5 records."""
+
+import importlib.util
+import json
+from pathlib import Path
+import subprocess
+import sys
+
+import h5py
+import numpy as np
+import pytest
+
+
+_SCRIPT = Path(__file__).parents[2] / 'tools/dev/compare_tracks.py'
+_SPEC = importlib.util.spec_from_file_location('compare_tracks', _SCRIPT)
+_TOOL = importlib.util.module_from_spec(_SPEC)
+_SPEC.loader.exec_module(_TOOL)
+
+_FLOAT_FIELDS = (
+    'r.x', 'r.y', 'r.z', 'u.x', 'u.y', 'u.z', 'E', 'time', 'wgt',
+)
+_INT_FIELDS = ('cell_id', 'cell_instance', 'material_id')
+
+
+def _dtype(byteorder='<', padded=False, reversed_fields=False):
+    xyz = [(axis, byteorder + 'f8') for axis in 'xyz']
+    fields = [('r', xyz), ('u', xyz)]
+    fields.extend((name, byteorder + 'f8') for name in ('E', 'time', 'wgt'))
+    fields.extend((name, byteorder + 'i4') for name in _INT_FIELDS)
+    if reversed_fields:
+        fields.reverse()
+    return np.dtype(fields, align=padded)
+
+
+def _field(states, field):
+    for name in field.split('.'):
+        states = states[name]
+    return states
+
+
+def _states(size=3, **dtype_kwargs):
+    states = np.zeros(size, dtype=_dtype(**dtype_kwargs))
+    for index, field in enumerate(_FLOAT_FIELDS):
+        _field(states, field)[:] = np.arange(size) + index + 0.5
+    for index, field in enumerate(_INT_FIELDS):
+        states[field] = np.arange(size) + index + 1
+    return states
+
+
+def _write(path, histories=None):
+    if histories is None:
+        histories = [('track_1_1_2', [(_states(), 2112)])]
+    with h5py.File(path, 'w') as fh:
+        fh.attrs['filetype'] = np.bytes_('track')
+        fh.attrs['version'] = np.array([3, 1], dtype='i4')
+        for name, particles in histories:
+            states = np.concatenate(
+                [state for state, _ in particles], dtype=particles[0][0].dtype)
+            dset = fh.create_dataset(name, data=states)
+            dset.attrs['n_particles'] = len(particles)
+            dset.attrs['particles'] = [pdg for _, pdg in particles]
+            dset.attrs['offsets'] = np.concatenate([
+                [0], np.cumsum([len(state) for state, _ in particles]),
+            ])
+    return path
+
+
+@pytest.fixture
+def track_pair(tmp_path):
+    return _write(tmp_path / 'left.h5'), _write(tmp_path / 'right.h5')
+
+
+def _compare_states(tmp_path, left, right):
+    a = _write(tmp_path / 'a.h5', [('track_1_1_2', [(left, 2112)])])
+    b = _write(tmp_path / 'b.h5', [('track_1_1_2', [(right, 2112)])])
+    return _TOOL.compare_tracks(a, b)
+
+
+def test_match(track_pair):
+    report = _TOOL.compare_tracks(*track_pair)
+    assert report == {
+        'status': 'match',
+        'left': {'histories': 1, 'particles': 1, 'states': 3},
+        'right': {'histories': 1, 'particles': 1, 'states': 3},
+    }
+
+
+@pytest.mark.parametrize('legacy,pdg', [(0, 2112), (1, 22), (2, 11), (3, -11)])
+def test_legacy_and_current_particle_codes_match(tmp_path, legacy, pdg):
+    a = _write(tmp_path / 'a.h5', [('track_1_1_2', [(_states(), legacy)])])
+    b = _write(tmp_path / 'b.h5', [('track_1_1_2', [(_states(), pdg)])])
+    with h5py.File(a, 'r+') as fh:
+        fh.attrs['version'] = [3, 0]
+    assert _TOOL.compare_tracks(a, b)['status'] == 'match'
+    assert _TOOL.compare_tracks(b, a)['status'] == 'match'
+
+
+def test_legacy_particle_type_difference(tmp_path):
+    paths = []
+    for legacy in (0, 1):
+        path = _write(tmp_path / f'{legacy}.h5', [
+            ('track_1_1_2', [(_states(), legacy)]),
+        ])
+        with h5py.File(path, 'r+') as fh:
+            fh.attrs['version'] = [3, 0]
+        paths.append(path)
+    report = _TOOL.compare_tracks(*paths)
+    assert report['status'] == 'mismatch'
+    assert report['difference']['kind'] == 'particle_type'
+    assert report['difference']['left'] == 2112
+    assert report['difference']['right'] == 22
+
+
+def test_legacy_file_rejects_mixed_particle_conventions(track_pair):
+    _write(track_pair[1], [
+        ('track_1_1_2', [(_states(), 0), (_states(), 22)]),
+    ])
+    with h5py.File(track_pair[1], 'r+') as fh:
+        fh.attrs['version'] = [3, 0]
+    report = _TOOL.compare_tracks(*track_pair)
+    assert report['status'] == 'invalid_input'
+    assert report['side'] == 'right'
+
+
+@pytest.mark.parametrize('field', _FLOAT_FIELDS + _INT_FIELDS)
+def test_each_state_field(tmp_path, field):
+    left, right = _states(), _states()
+    _field(right, field)[1] += 1
+    report = _compare_states(tmp_path, left, right)
+    assert report['status'] == 'mismatch'
+    difference = report['difference']
+    assert difference['history'] == [1, 1, 2]
+    assert difference['particle_index'] == 0
+    assert difference['particle_type'] == 2112
+    assert difference['kind'] == 'state'
+    assert difference['state_index'] == 1
+    assert difference['field'] == field
+    if field in _FLOAT_FIELDS:
+        assert difference['left']['bits'].startswith('0x')
+        assert difference['right']['bits'].startswith('0x')
+        assert difference['ulp_distance'] > 0
+    else:
+        assert difference['left']['value'] + 1 == difference['right']['value']
+
+
+def test_earliest_state_before_field(tmp_path):
+    left, right = _states(), _states()
+    right['r']['x'][2] += 1
+    right['material_id'][0] += 1
+    difference = _compare_states(tmp_path, left, right)['difference']
+    assert difference['state_index'] == 0
+    assert difference['field'] == 'material_id'
+
+
+def test_same_state_field_order(tmp_path):
+    left, right = _states(), _states()
+    right['u']['y'][1] += 1
+    right['u']['z'][1] += 1
+    right['E'][1] += 1
+    difference = _compare_states(tmp_path, left, right)['difference']
+    assert (difference['state_index'], difference['field']) == (1, 'u.y')
+
+
+@pytest.mark.parametrize('left_bits,right_bits,ulp', [
+    (0x0000000000000000, 0x8000000000000000, 0),
+    (0x0000000000000000, 0x0000000000000001, 1),
+    (0x8000000000000001, 0x0000000000000001, 2),
+    (0x0010000000000000, 0x000fffffffffffff, 1),
+    (0x3ff0000000000000, 0x3ff0000000000001, 1),
+    (0xbff0000000000000, 0xbff0000000000001, 1),
+    (0xbff0000000000000, 0x3ff0000000000000,
+     2 * 0x3ff0000000000000),
+    (0x7fefffffffffffff, 0xffefffffffffffff,
+     2 * 0x7fefffffffffffff),
+    (0x7ff0000000000000, 0xfff0000000000000, None),
+    (0x7fefffffffffffff, 0x7ff0000000000000, None),
+    (0x7ff8000000000001, 0x7ff8000000000002, None),
+    (0x7ff0000000000001, 0x7ff8000000000001, None),
+    (0x7ff8000000000001, 0xfff8000000000001, None),
+])
+def test_exact_float_bits(tmp_path, left_bits, right_bits, ulp):
+    left, right = _states(), _states()
+    left['E'].view('<u8')[1] = left_bits
+    right['E'].view('<u8')[1] = right_bits
+    report = _compare_states(tmp_path, left, right)
+    assert report['status'] == 'mismatch'
+    difference = report['difference']
+    assert difference['field'] == 'E'
+    assert difference['left']['bits'] == f'0x{left_bits:016x}'
+    assert difference['right']['bits'] == f'0x{right_bits:016x}'
+    assert difference['ulp_distance'] == ulp
+    # Nonfinite values must remain representable in standards-compliant JSON.
+    assert json.loads(json.dumps(report, allow_nan=False)) == report
+
+
+@pytest.mark.parametrize('bits', [
+    0x8000000000000000, 0x7ff0000000000000,
+    0xfff0000000000000, 0x7ff8000000000001, 0x7ff0000000000001,
+])
+def test_identical_nonfinite_and_signed_zero(tmp_path, bits):
+    left, right = _states(), _states()
+    left['E'].view('<u8')[1] = bits
+    right['E'].view('<u8')[1] = bits
+    assert _compare_states(tmp_path, left, right)['status'] == 'match'
+
+
+def test_storage_byteorder_padding_and_member_order(tmp_path):
+    left = _states()
+    right = _states(byteorder='>', padded=True, reversed_fields=True)
+    assert left.dtype != right.dtype
+    assert left.dtype.itemsize != right.dtype.itemsize
+    report = _compare_states(tmp_path, left, right)
+    assert report['status'] == 'match'
+    with h5py.File(tmp_path / 'b.h5', 'r') as fh:
+        dtype = fh['track_1_1_2'].dtype
+        assert dtype.names != left.dtype.names
+        assert dtype['E'].byteorder == '>'
+        assert dtype.itemsize != left.dtype.itemsize
+
+
+def test_big_endian_nan_payload(tmp_path):
+    left = _states()
+    right = _states(byteorder='>')
+    left['E'].view('<u8')[1] = 0x7ff8000000000001
+    right['E'].view('>u8')[1] = 0x7ff8000000000001
+    assert _compare_states(tmp_path, left, right)['status'] == 'match'
+    right['E'].view('>u8')[1] = 0x7ff8000000000002
+    difference = _compare_states(tmp_path, left, right)['difference']
+    assert difference['right']['bits'] == '0x7ff8000000000002'
+
+
+def test_numeric_history_order(tmp_path):
+    left, right = _states(), _states()
+    right['E'][0] += 1
+    histories = [
+        ('track_1_1_10', [(left, 2112)]),
+        ('track_1_1_2', [(left, 2112)]),
+    ]
+    a = _write(tmp_path / 'a.h5', histories)
+    b = _write(tmp_path / 'b.h5', [
+        (name, [(right, 2112)]) for name, _ in reversed(histories)
+    ])
+    assert _TOOL.compare_tracks(a, b)['difference']['history'] == [1, 1, 2]
+
+
+def test_particle_ordinal_before_state(tmp_path):
+    left, right = _states(), _states()
+    right['E'][2] += 1
+    a = _write(tmp_path / 'a.h5', [
+        ('track_1_1_2', [(left, 22), (left, 11)]),
+    ])
+    b = _write(tmp_path / 'b.h5', [
+        ('track_1_1_2', [(right, 22), (right, 11)]),
+    ])
+    difference = _TOOL.compare_tracks(a, b)['difference']
+    assert difference['particle_index'] == 0
+    assert difference['particle_type'] == 22
+    assert difference['state_index'] == 2
+
+
+@pytest.mark.parametrize('missing_from', ['left', 'right'])
+def test_missing_history(tmp_path, missing_from):
+    a = _write(tmp_path / 'a.h5')
+    b = _write(tmp_path / 'b.h5', [
+        ('track_1_1_2', [(_states(), 2112)]),
+        ('track_1_1_10', [(_states(), 2112)]),
+    ])
+    paths = (a, b) if missing_from == 'left' else (b, a)
+    report = _TOOL.compare_tracks(*paths)
+    assert report['status'] == 'mismatch'
+    assert report['difference'] == {
+        'history': [1, 1, 10], 'kind': 'missing_history',
+        'missing_from': missing_from,
+    }
+
+
+@pytest.mark.parametrize('kind', ['particle_type', 'state_count',
+                                 'particle_count'])
+def test_structural_difference(tmp_path, kind):
+    left = [(_states(), 2112), (_states(), 22)]
+    right = [(_states(), 2112), (_states(), 22)]
+    if kind == 'particle_type':
+        right[1] = (_states(), -11)
+    elif kind == 'state_count':
+        right[1] = (_states(2), 22)
+    else:
+        right.pop()
+    a = _write(tmp_path / 'a.h5', [('track_1_1_2', left)])
+    b = _write(tmp_path / 'b.h5', [('track_1_1_2', right)])
+    report = _TOOL.compare_tracks(a, b)
+    assert report['status'] == 'mismatch'
+    difference = report['difference']
+    assert difference['kind'] == kind
+    assert difference['particle_index'] == 1
+    if kind == 'state_count':
+        assert difference['state_index'] == 2
+        assert (difference['left'], difference['right']) == (3, 2)
+    elif kind == 'particle_count':
+        assert (difference['left'], difference['right']) == (2, 1)
+    else:
+        assert (difference['left'], difference['right']) == (22, -11)
+        assert 'particle_type' not in difference
+
+
+def test_state_difference_before_count_difference(tmp_path):
+    left, right = _states(), _states(2)
+    right['time'][0] += 1
+    difference = _compare_states(tmp_path, left, right)['difference']
+    assert difference['kind'] == 'state'
+    assert difference['state_index'] == 0
+    assert difference['field'] == 'time'
+
+
+@pytest.mark.parametrize('attribute,value', [
+    ('offsets', [-1, 3]),
+    ('offsets', [1, 3]),
+    ('offsets', [0, 2]),
+    ('offsets', [0, 4]),
+    ('offsets', [0]),
+    ('offsets', [0, 2, 3]),
+    ('offsets', [0.0, 3.0]),
+    ('offsets', [[0, 3]]),
+    ('offsets', np.array([0, 2**63], dtype='u8')),
+    ('particles', [0]),
+    ('particles', [1]),
+    ('particles', [2]),
+    ('particles', [3]),
+    ('particles', [2112.0]),
+    ('particles', [[2112]]),
+    ('particles', []),
+    ('n_particles', -1),
+    ('n_particles', 0),
+    ('n_particles', 2),
+    ('n_particles', [1]),
+    ('n_particles', 1.0),
+])
+def test_invalid_particle_metadata(track_pair, attribute, value):
+    with h5py.File(track_pair[1], 'r+') as fh:
+        fh['track_1_1_2'].attrs[attribute] = value
+    report = _TOOL.compare_tracks(*track_pair)
+    assert report['status'] == 'invalid_input'
+    assert report['side'] == 'right'
+
+
+def test_nonmonotone_unsigned_offsets(track_pair):
+    with h5py.File(track_pair[1], 'r+') as fh:
+        attrs = fh['track_1_1_2'].attrs
+        attrs['n_particles'] = 3
+        attrs['particles'] = [2112, 22, 11]
+        attrs['offsets'] = np.array([0, 2, 1, 3], dtype='u8')
+    assert _TOOL.compare_tracks(*track_pair)['status'] == 'invalid_input'
+
+
+@pytest.mark.parametrize('attribute', ['n_particles', 'particles', 'offsets'])
+def test_missing_particle_metadata(track_pair, attribute):
+    with h5py.File(track_pair[0], 'r+') as fh:
+        del fh['track_1_1_2'].attrs[attribute]
+    report = _TOOL.compare_tracks(*track_pair)
+    assert report['status'] == 'invalid_input'
+    assert report['side'] == 'left'
+
+
+@pytest.mark.parametrize('attribute,value', [
+    ('version', [2, 0]),
+    ('version', [4, 0]),
+    ('version', [3, 2]),
+    ('version', [3, -1]),
+    ('version', [3]),
+    ('version', [3, 1, 0]),
+    ('version', [3.0, 1.0]),
+    ('version', []),
+    ('version', [[3, 1]]),
+    ('version', 3),
+    ('filetype', np.bytes_('statepoint')),
+])
+def test_invalid_file_metadata(track_pair, attribute, value):
+    with h5py.File(track_pair[1], 'r+') as fh:
+        fh.attrs[attribute] = value
+    assert _TOOL.compare_tracks(*track_pair)['status'] == 'invalid_input'
+
+
+@pytest.mark.parametrize('attribute', ['version', 'filetype'])
+def test_missing_file_metadata(track_pair, attribute):
+    with h5py.File(track_pair[1], 'r+') as fh:
+        del fh.attrs[attribute]
+    assert _TOOL.compare_tracks(*track_pair)['status'] == 'invalid_input'
+
+
+@pytest.mark.parametrize('malformation', [
+    'group', 'scalar', 'matrix', 'wrong_float_width', 'wrong_integer_sign',
+    'missing_field', 'extra_field', 'coordinate_array', 'bad_name',
+    'duplicate_identifier', 'soft_link', 'external_link',
+])
+def test_invalid_dataset_schema(track_pair, malformation):
+    with h5py.File(track_pair[1], 'r+') as fh:
+        original = fh['track_1_1_2']
+        states = original[()]
+        attrs = dict(original.attrs)
+        if malformation == 'bad_name':
+            fh.move('track_1_1_2', 'track_one_1_2')
+        elif malformation == 'duplicate_identifier':
+            fh['track_01_1_2'] = original
+        elif malformation == 'soft_link':
+            fh['track_1_1_3'] = h5py.SoftLink('/track_1_1_2')
+        elif malformation == 'external_link':
+            fh['track_1_1_3'] = h5py.ExternalLink(
+                str(track_pair[0]), '/track_1_1_2')
+        else:
+            del fh['track_1_1_2']
+            if malformation == 'group':
+                fh.create_group('track_1_1_2')
+            elif malformation == 'scalar':
+                fh.create_dataset('track_1_1_2', data=states[0])
+            elif malformation == 'matrix':
+                fh.create_dataset('track_1_1_2', data=states.reshape(1, 3))
+            else:
+                fields = states.dtype.descr
+                if malformation == 'wrong_float_width':
+                    fields = [(n, '<f4' if n == 'E' else t) for n, t in fields]
+                elif malformation == 'wrong_integer_sign':
+                    fields = [(n, '<u4' if n == 'cell_id' else t)
+                              for n, t in fields]
+                elif malformation == 'missing_field':
+                    fields = [(n, t) for n, t in fields if n != 'time']
+                elif malformation == 'extra_field':
+                    fields.append(('extra', '<f8'))
+                elif malformation == 'coordinate_array':
+                    fields = [(n, ('<f8', (3,)) if n == 'r' else t)
+                              for n, t in fields]
+                dset = fh.create_dataset(
+                    'track_1_1_2', data=np.zeros(3, dtype=fields))
+                dset.attrs.update(attrs)
+    assert _TOOL.compare_tracks(*track_pair)['status'] == 'invalid_input'
+
+
+@pytest.mark.parametrize('empty_side', ['left', 'right', 'both'])
+def test_empty_file_is_not_agreement(track_pair, empty_side):
+    for index, side in enumerate(('left', 'right')):
+        if empty_side in (side, 'both'):
+            _write(track_pair[index], [])
+    report = _TOOL.compare_tracks(*track_pair)
+    if empty_side == 'both':
+        assert report['status'] == 'invalid_input'
+    else:
+        assert report['status'] == 'mismatch'
+        assert report['difference']['kind'] == 'missing_history'
+        assert report['difference']['missing_from'] == empty_side
+
+
+def test_late_malformed_history_takes_precedence(tmp_path):
+    left, right = _states(), _states()
+    right['E'][0] += 1
+    a = _write(tmp_path / 'a.h5')
+    b = _write(tmp_path / 'b.h5', [
+        ('track_1_1_2', [(right, 2112)]),
+        ('track_1_1_10', [(left, 2112)]),
+    ])
+    with h5py.File(b, 'r+') as fh:
+        fh['track_1_1_10'].attrs['offsets'] = [-1, 3]
+    report = _TOOL.compare_tracks(a, b)
+    assert report['status'] == 'invalid_input'
+    assert report['side'] == 'right'
+
+
+def test_empty_particle_states_are_not_agreement(tmp_path):
+    report = _compare_states(tmp_path, _states(0), _states(0))
+    assert report['status'] == 'invalid_input'
+
+
+@pytest.mark.parametrize('storage', ['external_link', 'external', 'virtual'])
+def test_external_storage_rejected_before_read(track_pair, tmp_path, storage):
+    with h5py.File(track_pair[1], 'r+') as fh:
+        attrs = dict(fh['track_1_1_2'].attrs)
+        del fh['track_1_1_2']
+        absent = str(tmp_path / 'not-present.h5')
+        if storage == 'external_link':
+            fh['track_1_1_2'] = h5py.ExternalLink(absent, '/records')
+        elif storage == 'external':
+            dset = fh.create_dataset(
+                'track_1_1_2', shape=(3,), dtype=_dtype(),
+                external=[(absent, 0, h5py.h5f.UNLIMITED)],
+            )
+            dset.attrs.update(attrs)
+        else:
+            layout = h5py.VirtualLayout(shape=(3,), dtype=_dtype())
+            layout[:] = h5py.VirtualSource(absent, 'records', shape=(3,))
+            dset = fh.create_virtual_dataset('track_1_1_2', layout)
+            dset.attrs.update(attrs)
+    report = _TOOL.compare_tracks(*track_pair)
+    assert report['status'] == 'invalid_input'
+    expected = 'linked track data' if storage == 'external_link' else (
+        'self-contained track data')
+    assert expected in report['error']
+
+
+def test_comparison_does_not_modify_inputs(track_pair):
+    before = [path.read_bytes() for path in track_pair]
+    assert _TOOL.compare_tracks(*track_pair)['status'] == 'match'
+    assert [path.read_bytes() for path in track_pair] == before
+
+
+def test_missing_file(track_pair, tmp_path):
+    report = _TOOL.compare_tracks(track_pair[0], tmp_path / 'missing.h5')
+    assert report['status'] == 'invalid_input'
+    assert report['side'] == 'right'
+
+
+@pytest.mark.parametrize('status,code', [
+    ('match', 0), ('mismatch', 1), ('invalid_input', 2),
+])
+def test_cli_exit_status(track_pair, tmp_path, status, code):
+    left, right = track_pair
+    if status == 'mismatch':
+        states = _states()
+        states['E'][0] += 1
+        _write(right, [('track_1_1_2', [(states, 2112)])])
+    elif status == 'invalid_input':
+        right = tmp_path / 'missing.h5'
+    result = subprocess.run(
+        [sys.executable, str(_SCRIPT), str(left), str(right)],
+        text=True, capture_output=True, check=False,
+    )
+    assert result.returncode == code, result.stderr
+    assert json.loads(result.stdout)['status'] == status
+    assert 'Traceback' not in result.stderr
diff --git a/tools/dev/compare_tracks.py b/tools/dev/compare_tracks.py
new file mode 100644
index 00000000000..2e360fc6f1f
--- /dev/null
+++ b/tools/dev/compare_tracks.py
@@ -0,0 +1,231 @@
+#!/usr/bin/env python3
+"""Report the first exact difference between two OpenMC track files."""
+
+import argparse
+import json
+from pathlib import Path
+import re
+
+import h5py
+import numpy as np
+import openmc
+
+
+# Logical order, independent of compound-member offsets and HDF5 byte order.
+_FIELDS = (
+    'r.x', 'r.y', 'r.z', 'u.x', 'u.y', 'u.z', 'E', 'time', 'wgt',
+    'cell_id', 'cell_instance', 'material_id',
+)
+_FLOAT_FIELDS = _FIELDS[:9]
+_SIGN = 1 << 63
+_EXPONENT = 0x7ff0000000000000
+
+
+def _field(array, name):
+    for component in name.split('.'):
+        array = array[component]
+    return array
+
+
+def _validate_dtype(dtype, fields):
+    """Validate typed members, not padding, offsets, or byte order."""
+    if dtype.names is None or set(dtype.names) != set(fields):
+        raise ValueError(f'Expected state fields {tuple(fields)}')
+    for name, expected in fields.items():
+        member = dtype.fields[name][0]
+        if isinstance(expected, dict):
+            _validate_dtype(member, expected)
+        elif member.subdtype or (member.kind, member.itemsize) != expected:
+            raise ValueError(f'Invalid state dtype for {name}: {member}')
+
+
+def _load_tracks(path):
+    """Validate structural metadata before the Tracks reader slices records."""
+    xyz = dict.fromkeys(('x', 'y', 'z'), ('f', 8))
+    fields = {'r': xyz, 'u': xyz, 'E': ('f', 8), 'time': ('f', 8),
+              'wgt': ('f', 8), 'cell_id': ('i', 4),
+              'cell_instance': ('i', 4), 'material_id': ('i', 4)}
+    identifiers = set()
+    with h5py.File(path, 'r') as fh:
+        version = np.asarray(fh.attrs['version'])
+        if (version.ndim != 1 or version.size != 2
+                or version.dtype.kind not in 'iu'
+                or tuple(version) not in ((3, 0), (3, 1))):
+            raise ValueError('Expected track format version 3.0 or 3.1')
+        openmc.checkvalue.check_filetype_version(fh, 'track', 3)
+        for name in fh:
+            match = re.fullmatch(r'track_([0-9]+)_([0-9]+)_([0-9]+)', name)
+            if match is None:
+                raise ValueError(f'Invalid history name: {name}')
+            identifier = tuple(map(int, match.groups()))
+            if identifier in identifiers:
+                raise ValueError(f'Duplicate history identifier: {identifier}')
+            identifiers.add(identifier)
+            if not isinstance(fh.get(name, getlink=True), h5py.HardLink):
+                raise ValueError(f'{name}: linked track data are not supported')
+            dset = fh[name]
+            if not isinstance(dset, h5py.Dataset) or dset.ndim != 1:
+                raise ValueError(f'{name}: expected a one-dimensional dataset')
+            if dset.is_virtual or dset.external:
+                raise ValueError(f'{name}: expected self-contained track data')
+            _validate_dtype(dset.dtype, fields)
+            particles = np.asarray(dset.attrs['particles'])
+            offsets = np.asarray(dset.attrs['offsets'])
+            count = np.asarray(dset.attrs['n_particles'])
+            if (count.ndim != 0 or count.dtype.kind not in 'iu'
+                    or particles.ndim != 1 or particles.dtype.kind not in 'iu'
+                    or int(count) != particles.size or particles.size == 0):
+                raise ValueError(f'{name}: inconsistent particle count')
+            # Native tracks include at least their final state. Comparing
+            # adjacent entries also avoids unsigned subtraction overflow.
+            if (offsets.ndim != 1 or offsets.dtype.kind not in 'iu'
+                    or offsets.size != particles.size + 1
+                    or offsets[0] != 0 or offsets[-1] != len(dset)
+                    or np.any(offsets[1:] <= offsets[:-1])):
+                raise ValueError(f'{name}: invalid particle offsets')
+            # Track 3.0 used legacy indices; 3.1 uses PDG numbers. The Tracks
+            # reader normalizes the former to PDG. Validate the encoding first
+            # so that this normalization cannot hide malformed 3.1 metadata.
+            legacy = np.isin(particles, [0, 1, 2, 3])
+            if version[1] == 0 and not np.all(legacy):
+                raise ValueError(f'{name}: expected legacy particle indices')
+            if version[1] == 1 and np.any(legacy):
+                raise ValueError(f'{name}: expected PDG particle numbers')
+    return openmc.Tracks(path)
+
+
+def _bits(values):
+    """View binary64 bits without floating arithmetic or NaN canonicalization."""
+    return values.view(np.dtype('u8').newbyteorder(values.dtype.byteorder))
+
+
+def _value(value, floating):
+    if not floating:
+        return {'value': int(value)}
+    bits = int(_bits(np.asarray(value)))
+    return {'value': repr(float(value)), 'bits': f'0x{bits:016x}'}
+
+
+def _ulp_distance(left, right):
+    """Count finite representable steps, with the two signed zeros coalesced."""
+    if left & _EXPONENT == _EXPONENT or right & _EXPONENT == _EXPONENT:
+        return None
+    left_rank = -(left & (_SIGN - 1)) if left & _SIGN else left
+    right_rank = -(right & (_SIGN - 1)) if right & _SIGN else right
+    return abs(left_rank - right_rank)
+
+
+def _first_state_difference(left, right):
+    """Find the earliest state, breaking ties in documented field order."""
+    count = min(len(left), len(right))
+    first = None
+    for order, name in enumerate(_FIELDS):
+        a, b = _field(left[:count], name), _field(right[:count], name)
+        floating = name in _FLOAT_FIELDS
+        aa, bb = (_bits(a), _bits(b)) if floating else (a, b)
+        different = np.flatnonzero(aa != bb)
+        if different.size:
+            index = int(different[0])
+            if first is None or (index, order) < first[:2]:
+                detail = {'kind': 'state', 'state_index': index, 'field': name,
+                          'left': _value(a[index], floating),
+                          'right': _value(b[index], floating)}
+                if floating:
+                    detail['ulp_distance'] = _ulp_distance(
+                        int(aa[index]), int(bb[index]))
+                first = (index, order, detail)
+    return None if first is None else first[2]
+
+
+def _summary(tracks):
+    return {'histories': len(tracks),
+            'particles': sum(len(t) for t in tracks),
+            'states': sum(len(p.states) for t in tracks for p in t)}
+
+
+def compare_tracks(left, right):
+    """Compare recorded histories, returning a JSON-serializable report.
+
+    Parameters
+    ----------
+    left, right : str or pathlib.Path
+        Self-contained version-3.0 or 3.1 track files. Both must be closed by the
+        producing simulations before comparison.
+
+    Returns
+    -------
+    dict
+        ``status`` is ``match``, ``mismatch``, or ``invalid_input``. On mismatch,
+        ``difference`` identifies the first difference in history identifier,
+        particle-track sequence, state index, then field order. Secondary
+        indices are recorded sequence positions, not stable genealogy IDs.
+
+    """
+    inputs = []
+    for side, path in (('left', left), ('right', right)):
+        try:
+            inputs.append(_load_tracks(path))
+        except (OSError, ValueError, TypeError, KeyError, IndexError,
+                OverflowError) as exc:
+            return {'status': 'invalid_input', 'side': side, 'error': str(exc)}
+
+    a_tracks, b_tracks = inputs
+    report = {'status': 'match', 'left': _summary(a_tracks),
+              'right': _summary(b_tracks)}
+    if not a_tracks and not b_tracks:
+        # Empty files cannot establish agreement of any recorded history.
+        return {'status': 'invalid_input',
+                'error': 'Neither input contains a recorded history'}
+    a_by_id = {t.identifier: t for t in a_tracks}
+    b_by_id = {t.identifier: t for t in b_tracks}
+
+    for identifier in sorted(a_by_id.keys() | b_by_id.keys()):
+        a, b = a_by_id.get(identifier), b_by_id.get(identifier)
+        location = {'history': list(identifier)}
+        difference = None
+        if a is None or b is None:
+            difference = {'kind': 'missing_history',
+                          'missing_from': 'left' if a is None else 'right'}
+        else:
+            for index, (pa, pb) in enumerate(zip(a, b)):
+                location = {'history': list(identifier),
+                            'particle_index': index}
+                if pa.particle != pb.particle:
+                    difference = {'kind': 'particle_type',
+                                  'left': int(pa.particle),
+                                  'right': int(pb.particle)}
+                else:
+                    location['particle_type'] = int(pa.particle)
+                    difference = _first_state_difference(pa.states, pb.states)
+                    if difference is None and len(pa.states) != len(pb.states):
+                        difference = {'kind': 'state_count',
+                                      'state_index': min(len(pa.states),
+                                                         len(pb.states)),
+                                      'left': len(pa.states),
+                                      'right': len(pb.states)}
+                if difference is not None:
+                    break
+            if difference is None and len(a) != len(b):
+                location = {'history': list(identifier),
+                            'particle_index': min(len(a), len(b))}
+                difference = {'kind': 'particle_count',
+                              'left': len(a), 'right': len(b)}
+        if difference is not None:
+            report['status'] = 'mismatch'
+            report['difference'] = {**location, **difference}
+            return report
+    return report
+
+
+def main(argv=None):
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument('left', type=Path)
+    parser.add_argument('right', type=Path)
+    args = parser.parse_args(argv)
+    report = compare_tracks(args.left, args.right)
+    print(json.dumps(report, indent=2, allow_nan=False))
+    return {'match': 0, 'mismatch': 1, 'invalid_input': 2}[report['status']]
+
+
+if __name__ == '__main__':
+    raise SystemExit(main())