From 8cf6b1a5c657a308644d76ceecd6ec4233c3fef8 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 15 Aug 2026 16:18:41 -0700 Subject: [PATCH 1/4] Forward all supported params to the Fortran run --- pamica/tests/test_fortran_param_forwarding.py | 107 ++++++++++++++++++ validate_implementations.py | 103 ++++++++++++++--- 2 files changed, 193 insertions(+), 17 deletions(-) create mode 100644 pamica/tests/test_fortran_param_forwarding.py diff --git a/pamica/tests/test_fortran_param_forwarding.py b/pamica/tests/test_fortran_param_forwarding.py new file mode 100644 index 00000000..de8951ff --- /dev/null +++ b/pamica/tests/test_fortran_param_forwarding.py @@ -0,0 +1,107 @@ +"""Parity-harness parameter forwarding (issue #228). + +Both arms of a parity run must be configured identically. The writer previously +rewrote six hardcoded keys and left everything else at the template's value while +the Python side honored it, which silently makes the comparison uncontrolled. +""" + +import sys +from pathlib import Path + +import pytest + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT)) + +from validate_implementations import ( # noqa: E402 + fortran_accepted_keys, + write_fortran_param_file, +) + +TEMPLATE = ROOT / "pamica" / "sample_data" / "input.param" +SOURCE = ROOT / "pamica" / "amica15.f90" + +pytestmark = pytest.mark.skipif( + not TEMPLATE.exists() or not SOURCE.exists(), + reason="reference param template or Fortran source missing", +) + + +def _written(tmp_path, params, overrides=None): + dest = tmp_path / "input.param" + write_fortran_param_file( + TEMPLATE.read_text().splitlines(keepends=True), dest, params, overrides + ) + lines = [ln for ln in dest.read_text().splitlines() if ln.split()] + out = dict(ln.split(None, 1) for ln in lines) + assert len(out) == len(lines), "a key was written more than once" + return {k: v.strip() for k, v in out.items()} + + +def test_accepted_keys_come_from_the_reference_source(): + keys = fortran_accepted_keys(SOURCE) + assert keys is not None + # Spot-check across the parser rather than pinning the whole set, which + # would just restate the source. + for key in ("do_newton", "block_size", "do_approx_sphere", "max_decs", "lrate"): + assert key in keys, f"{key} should be parsed out of amica15.f90" + + +def test_missing_source_disables_filtering(): + assert fortran_accepted_keys(Path("does-not-exist.f90")) is None + + +def test_do_newton_is_forwarded(tmp_path): + """The regression that motivated #228: this key was silently dropped, so + ``params['do_newton'] = False`` configured Python without Newton and left the + reference running with it.""" + assert _written(tmp_path, {"do_newton": False})["do_newton"] == "0" + assert _written(tmp_path, {"do_newton": True})["do_newton"] == "1" + + +def test_block_size_is_forwarded(tmp_path): + """Needed since the default moved to 8192 (#229) while the template says 512.""" + assert _written(tmp_path, {"block_size": 8192})["block_size"] == "8192" + + +def test_aliased_names_are_translated(tmp_path): + got = _written(tmp_path, {"num_mix": 4, "share_int": 50, "maxrej": 7}) + assert got["num_mix_comps"] == "4" + assert got["share_iter"] == "50" + assert got["numrej"] == "7" + + +def test_keys_absent_from_the_template_are_appended(tmp_path): + """The binary accepts more keywords than the shipped template lists.""" + assert "do_approx_sphere" not in TEMPLATE.read_text() + assert _written(tmp_path, {"do_approx_sphere": True})["do_approx_sphere"] == "1" + + +def test_overrides_beat_params(tmp_path): + """The harness's paths must win over whatever params.json carries.""" + got = _written( + tmp_path, + {"files": ["some/other/path.fdt"], "outdir": "./elsewhere/"}, + overrides={"files": "./eeglab_data.fdt", "outdir": "./fortran_output/"}, + ) + assert got["files"] == "./eeglab_data.fdt" + assert got["outdir"] == "./fortran_output/" + + +def test_untouched_template_values_survive(tmp_path): + got = _written(tmp_path, {"max_iter": 300}) + assert got["max_iter"] == "300" + assert got["data_dim"] == "32" # not in params, must keep the template's value + + +def test_unsupported_key_is_reported(tmp_path, capsys): + """A setting the reference cannot honor must never be dropped in silence.""" + _written(tmp_path, {"mineig_rel": 1e-12}) + assert "mineig_rel" in capsys.readouterr().out + + +def test_python_only_keys_are_not_reported(tmp_path, capsys): + """``device``/``seed`` configure our side only, so their absence is expected.""" + _written(tmp_path, {"device": "cpu", "seed": 42}) + out = capsys.readouterr().out + assert "device" not in out and "seed" not in out diff --git a/validate_implementations.py b/validate_implementations.py index b90288d4..f20142fe 100644 --- a/validate_implementations.py +++ b/validate_implementations.py @@ -9,6 +9,7 @@ import numpy as np import json import os +import re import subprocess import torch import inspect @@ -80,6 +81,85 @@ def load_sample_data() -> Tuple[np.ndarray, Dict]: return data, params +# params.json spellings that differ from the Fortran keyword of the same setting. +_FORTRAN_ALIASES = { + "num_mix": "num_mix_comps", + "share_int": "share_iter", + "maxrej": "numrej", +} + +# params.json keys that configure the Python side only, so having no Fortran +# keyword is expected rather than a dropped setting. +_PYTHON_ONLY_KEYS = {"device", "seed", "outdir", "files", "dtype"} + + +def fortran_accepted_keys( + source: Path = Path("pamica/amica15.f90"), +) -> Optional[set]: + """Keywords the reference binary's parameter parser accepts. + + Read from its own ``case('...')`` arms rather than hardcoded, so this cannot + drift from the binary being run. Returns ``None`` if the source is + unavailable, which callers treat as "forward everything and do not warn". + """ + if not source.exists(): + return None + return set(re.findall(r"^\s*case\('([^']+)'\)", source.read_text(), re.MULTILINE)) + + +def write_fortran_param_file( + template_lines, + dest: Path, + params: Dict, + overrides: Optional[Dict] = None, +) -> None: + """Write ``dest`` with every setting in ``params`` the binary understands. + + The Fortran and Python arms of a parity run must be configured identically. + Rewriting only a hardcoded handful of keys (as this did before issue #228) + left the rest at the template's values while the Python side honored them, + which silently turns a parity comparison into an uncontrolled one. + + Keys absent from the template are appended, since the binary accepts more + keywords than the shipped ``input.param`` lists. Keys it does not accept are + reported, so a setting can never be dropped in silence. + """ + accepted = fortran_accepted_keys() + wanted, unsupported = {}, [] + for key, value in params.items(): + fortran_key = _FORTRAN_ALIASES.get(key, key) + if accepted is not None and fortran_key not in accepted: + if key not in _PYTHON_ONLY_KEYS: + unsupported.append(key) + continue + # Fortran parses logicals as 0/1 integers. + wanted[fortran_key] = int(value) if isinstance(value, bool) else value + + # Last, so the harness's own paths win over whatever params.json carries. + wanted.update(overrides or {}) + + if unsupported: + print( + "WARNING: params with no Fortran keyword are left at the template's " + f"value for the reference run: {sorted(unsupported)}. The Python run " + "may honor them, so the comparison would not be controlled." + ) + + written, out = set(), [] + for line in template_lines: + key = line.split()[0] if line.split() else None + if key in wanted: + out.append(f"{key} {wanted[key]}\n") + written.add(key) + else: + out.append(line) + for key, value in wanted.items(): + if key not in written: + out.append(f"{key} {value}\n") + + dest.write_text("".join(out)) + + def run_fortran_amica( data: np.ndarray, params: Dict, @@ -120,23 +200,12 @@ def run_fortran_amica( with open(sample_param_file, "r") as f: param_lines = f.readlines() - # Update parameter file with our settings - with open(working_param_file, "w") as f: - for line in param_lines: - if line.startswith("files"): - f.write("files ./eeglab_data.fdt\n") - elif line.startswith("outdir"): - f.write("outdir ./fortran_output/\n") - elif line.startswith("max_iter"): - f.write(f"max_iter {params.get('max_iter', 100)}\n") - elif line.startswith("lrate"): - f.write(f"lrate {params.get('lrate', 0.05)}\n") - elif line.startswith("pdftype"): - f.write(f"pdftype {params.get('pdftype', 0)}\n") - elif line.startswith("num_mix_comps"): - f.write(f"num_mix_comps {params.get('num_mix', 3)}\n") - else: - f.write(line) + write_fortran_param_file( + param_lines, + working_param_file, + params, + overrides={"files": "./eeglab_data.fdt", "outdir": "./fortran_output/"}, + ) # Create output directory fortran_output = fortran_dir / "fortran_output" From 9f11e7706039bb18b4bf8a9a8b863579f595b92d Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 15 Aug 2026 16:46:14 -0700 Subject: [PATCH 2/4] Carry seedable Fortran in the tracked source --- native/patch_sources.py | 30 ++++++++++++++++++------------ pamica/amica15.f90 | 24 +++++++++++++++++++++++- pamica/amica15_header.f90 | 3 +++ 3 files changed, 44 insertions(+), 13 deletions(-) diff --git a/native/patch_sources.py b/native/patch_sources.py index cf4cd731..a465b611 100755 --- a/native/patch_sources.py +++ b/native/patch_sources.py @@ -8,8 +8,11 @@ an anchor is absent (source drift) it fails loudly rather than building a subtly wrong binary. -Only the build copy is patched; the tracked ``amica15.f90`` stays the read-only -parity reference (== sccn/amica master until PR #54 merges). +The tracked ``amica15.f90`` now carries PR #54 already (issue #228): pamica +maintains an improved copy rather than mirroring sccn/amica master, because the +unpatched reference cannot be seeded and so cannot anchor a parity gate. These +patches are therefore normally no-ops, and are kept so the script still works on +an unpatched upstream source. """ import sys @@ -81,21 +84,24 @@ def main(): pin_seed = "--pin-seed" in sys.argv[1:] amica15, header = args[0], args[1] - new_seed = NEW_SEED + patch(header, HDR_ANCHOR, HDR_ADD, "INTEGER :: nseed", "header seedvec decl") + patch(amica15, OLD_SEED, NEW_SEED, "call random_seed(size = nseed)", "random_seed") + patch(amica15, OLD_BCAST, NEW_BCAST, "MPI_BCAST(use_seed", "seed BCAST") + patch(amica15, OLD_PARSE, NEW_PARSE, "case('seed')", "seed param parser") + if pin_seed: # Determinism for validate_shim.sh: drop the clock term (c1) from the # default (unseeded) branch so the shim and mpif90 builds seed identically - # and can be compared bit-for-bit. Scoped to that one formula. - if "seedvec(jj) = c1 +" not in NEW_SEED: + # and can be compared bit-for-bit. Applied to the file after patching, not + # to NEW_SEED, so it still takes effect when the source already carries + # PR #54 and the patch above was a no-op. + text = open(amica15).read() + if "seedvec(jj) = c1 +" not in text: sys.exit( - "ERROR: --pin-seed anchor 'seedvec(jj) = c1 +' not found in NEW_SEED" + f"ERROR: --pin-seed anchor 'seedvec(jj) = c1 +' not found in {amica15}" ) - new_seed = NEW_SEED.replace("seedvec(jj) = c1 +", "seedvec(jj) = 987654321 +") - - patch(header, HDR_ANCHOR, HDR_ADD, "INTEGER :: nseed", "header seedvec decl") - patch(amica15, OLD_SEED, new_seed, "call random_seed(size = nseed)", "random_seed") - patch(amica15, OLD_BCAST, NEW_BCAST, "MPI_BCAST(use_seed", "seed BCAST") - patch(amica15, OLD_PARSE, NEW_PARSE, "case('seed')", "seed param parser") + with open(amica15, "w") as f: + f.write(text.replace("seedvec(jj) = c1 +", "seedvec(jj) = 987654321 +", 1)) print( "patched: portable+reproducible random_seed, seed param, BCAST" f"{' (pinned)' if pin_seed else ''}" diff --git a/pamica/amica15.f90 b/pamica/amica15.f90 index 10243bfd..1f1f2dd7 100644 --- a/pamica/amica15.f90 +++ b/pamica/amica15.f90 @@ -214,10 +214,28 @@ program main call MPI_BCAST(minlrate,1,MPI_DOUBLE_PRECISION,0,seg_comm,ierr) call MPI_BCAST(lratefact,1,MPI_DOUBLE_PRECISION,0,seg_comm,ierr) call MPI_BCAST(rholratefact,1,MPI_DOUBLE_PRECISION,0,seg_comm,ierr) +call MPI_BCAST(input_seed,1,MPI_INTEGER,0,seg_comm,ierr) +call MPI_BCAST(use_seed,1,MPI_LOGICAL,0,seg_comm,ierr) !--- set up the random number generator for this node call system_clock(c1) -call random_seed(PUT = c1 * (myrank+1) * (seed+myrank+1)) +! Portable + reproducible RNG seeding (sccn/amica PR #54, supersedes #53): +! random_seed(PUT=...) needs an array of the compiler SIZE (query with SIZE=); +! the original size-2 array only compiled under ifort. With `seed` set in the +! param file, seed deterministically (reproducible); otherwise clock-seed. +call random_seed(size = nseed) +allocate(seedvec(nseed)) +if (use_seed) then + do jj = 1, nseed + seedvec(jj) = input_seed + 1009*myrank + 37*(jj-1) + end do +else + do jj = 1, nseed + seedvec(jj) = c1 + 1009*(myrank+1)*jj + 37*(jj-1) + end do +end if +call random_seed(PUT = seedvec) +deallocate(seedvec) !call DRANDINITIALIZE(1,1,(myrank+1)*(c1/tot_procs + 1),lseed,state,lstate,info) @@ -3675,6 +3693,10 @@ subroutine get_cmd_args read(tmparg,'(a)') outdirparam case('indir') read(tmparg,'(a)') indirparam + case('seed') + read(tmparg,'(i12)') input_seed + use_seed = .true. + print *, 'seed = ', input_seed; call flush(6) end select end do diff --git a/pamica/amica15_header.f90 b/pamica/amica15_header.f90 index 89e66613..046f7afe 100644 --- a/pamica/amica15_header.f90 +++ b/pamica/amica15_header.f90 @@ -99,6 +99,9 @@ integer :: sampstop, blknum, blocknum, fld1, fld2, num_blocks INTEGER :: i, j, k integer :: ii, jj, kk, c0, c1, c2 +INTEGER :: nseed, input_seed = 0 +LOGICAL :: use_seed = .false. +INTEGER, ALLOCATABLE :: seedvec(:) integer :: counts_per_sec=0, cnt1=0, cnt2=0 INTEGER :: iter, len, fh, info, lwork, lstate=16, state(16), lseed = 1 integer :: host_num, ip(4), name_len From 16534d87c93c24407cddf6486aef1ec2752bc931 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 15 Aug 2026 16:48:35 -0700 Subject: [PATCH 3/4] Seed and pin threads for the reference run --- pamica/tests/test_fortran_param_forwarding.py | 25 ++++++++++++ validate_implementations.py | 40 ++++++++++++++++--- 2 files changed, 59 insertions(+), 6 deletions(-) diff --git a/pamica/tests/test_fortran_param_forwarding.py b/pamica/tests/test_fortran_param_forwarding.py index de8951ff..dee375d5 100644 --- a/pamica/tests/test_fortran_param_forwarding.py +++ b/pamica/tests/test_fortran_param_forwarding.py @@ -105,3 +105,28 @@ def test_python_only_keys_are_not_reported(tmp_path, capsys): _written(tmp_path, {"device": "cpu", "seed": 42}) out = capsys.readouterr().out assert "device" not in out and "seed" not in out + + +def test_real_params_json_produces_a_parseable_file(tmp_path): + """Regression: the full params.json, not a hand-picked subset. + + `files` and `field_dim` are lists there, and writing Python's repr put + brackets in the file, which aborts the Fortran parser at read time. The + unit tests above all passed while the harness could not run. + """ + import json + + params = json.loads( + (ROOT / "pamica" / "sample_data" / "sample_params.json").read_text() + ) + got = _written(tmp_path, params, overrides={"seed": 42, "max_threads": 1}) + for key, value in got.items(): + assert "[" not in value and "]" not in value, f"{key} kept Python list syntax" + assert got["field_dim"] == "30504" + assert got["seed"] == "42" + assert got["max_threads"] == "1" + + +def test_list_values_are_space_separated(tmp_path): + got = _written(tmp_path, {"field_dim": [100, 200]}) + assert got["field_dim"] == "100 200" diff --git a/validate_implementations.py b/validate_implementations.py index f20142fe..28a2ddfd 100644 --- a/validate_implementations.py +++ b/validate_implementations.py @@ -90,7 +90,13 @@ def load_sample_data() -> Tuple[np.ndarray, Dict]: # params.json keys that configure the Python side only, so having no Fortran # keyword is expected rather than a dropped setting. -_PYTHON_ONLY_KEYS = {"device", "seed", "outdir", "files", "dtype"} +_PYTHON_ONLY_KEYS = {"device", "outdir", "files", "dtype"} + +# Keywords the tracked amica15.f90 carries but sccn/amica master does not, so +# they are absent when this parses an upstream source (issue #228). The binary's +# parser has no `case default`, so an unknown keyword is silently ignored -- +# passing `seed` to an unpatched legacy binary is harmless, it just has no effect. +_PAMICA_EXTRA_KEYS = {"seed"} def fortran_accepted_keys( @@ -104,7 +110,22 @@ def fortran_accepted_keys( """ if not source.exists(): return None - return set(re.findall(r"^\s*case\('([^']+)'\)", source.read_text(), re.MULTILINE)) + found = set(re.findall(r"^\s*case\('([^']+)'\)", source.read_text(), re.MULTILINE)) + return found | _PAMICA_EXTRA_KEYS + + +def _fortran_value(value) -> str: + """Render a params.json value the way the Fortran parser reads it. + + Logicals are 0/1 integers, and the per-file lists (``files``, ``field_dim``) + are whitespace-separated -- writing Python's ``repr`` for those puts brackets + in the file and the parser aborts. + """ + if isinstance(value, bool): + return str(int(value)) + if isinstance(value, (list, tuple)): + return " ".join(_fortran_value(v) for v in value) + return str(value) def write_fortran_param_file( @@ -132,11 +153,10 @@ def write_fortran_param_file( if key not in _PYTHON_ONLY_KEYS: unsupported.append(key) continue - # Fortran parses logicals as 0/1 integers. - wanted[fortran_key] = int(value) if isinstance(value, bool) else value + wanted[fortran_key] = _fortran_value(value) # Last, so the harness's own paths win over whatever params.json carries. - wanted.update(overrides or {}) + wanted.update({k: _fortran_value(v) for k, v in (overrides or {}).items()}) if unsupported: print( @@ -200,11 +220,19 @@ def run_fortran_amica( with open(sample_param_file, "r") as f: param_lines = f.readlines() + # seed and max_threads are pinned, not taken from params: an unseeded or + # multi-threaded reference run is not reproducible (issue #228), which makes + # any comparison against it a comparison with a random draw. write_fortran_param_file( param_lines, working_param_file, params, - overrides={"files": "./eeglab_data.fdt", "outdir": "./fortran_output/"}, + overrides={ + "files": "./eeglab_data.fdt", + "outdir": "./fortran_output/", + "seed": int(seed), + "max_threads": 1, + }, ) # Create output directory From 1565346bcd9bb846d14331b3b08c9eb6a07e9fd2 Mon Sep 17 00:00:00 2001 From: Seyed Yahya Shirazi Date: Sat, 15 Aug 2026 16:56:30 -0700 Subject: [PATCH 4/4] Default the reference run to the native engine --- pamica/tests/test_fortran_param_forwarding.py | 26 +++++++++++ validate_implementations.py | 46 +++++++++++++++++-- 2 files changed, 68 insertions(+), 4 deletions(-) diff --git a/pamica/tests/test_fortran_param_forwarding.py b/pamica/tests/test_fortran_param_forwarding.py index dee375d5..c5244aba 100644 --- a/pamica/tests/test_fortran_param_forwarding.py +++ b/pamica/tests/test_fortran_param_forwarding.py @@ -130,3 +130,29 @@ def test_real_params_json_produces_a_parseable_file(tmp_path): def test_list_values_are_space_separated(tmp_path): got = _written(tmp_path, {"field_dim": [100, 200]}) assert got["field_dim"] == "100 200" + + +def test_default_reference_binary_falls_back_loudly(capsys, monkeypatch): + """Without a native engine the harness must say the run is uncontrolled. + + The fallback binary cannot be seeded, so silently using it would report a + comparison against a random draw as if it were controlled (issue #228). + """ + from validate_implementations import LEGACY_BINARY, default_reference_binary + + monkeypatch.delenv("PAMICA_NATIVE_BINARY", raising=False) + got = default_reference_binary(download=False) + if got == LEGACY_BINARY: + out = capsys.readouterr().out + assert "cannot be seeded" in out and "#228" in out + else: # a cached native engine is present; that is the good path + assert got.exists() + + +def test_env_override_selects_the_native_engine(tmp_path, monkeypatch): + fake = tmp_path / "amica15_native" + fake.write_text("#!/bin/sh\n") + monkeypatch.setenv("PAMICA_NATIVE_BINARY", str(fake)) + from validate_implementations import default_reference_binary + + assert default_reference_binary(download=False) == fake.resolve() diff --git a/validate_implementations.py b/validate_implementations.py index 28a2ddfd..1127e8ee 100644 --- a/validate_implementations.py +++ b/validate_implementations.py @@ -180,6 +180,33 @@ def write_fortran_param_file( dest.write_text("".join(out)) +LEGACY_BINARY = Path("pamica/sample_data/amica15mac") + + +def default_reference_binary(*, download: bool = True) -> Path: + """The reference binary to compare against, preferring the native engine. + + The native engine is built from a source carrying the ``seed`` option + (sccn/amica PR #54), so its runs are reproducible; the bundled + ``amica15mac`` fixture predates it and re-randomizes its initialization on + every run, which makes any comparison against it a comparison with a random + draw (issue #228). Falls back to the fixture, loudly, when the native engine + cannot be resolved. + """ + try: + from pamica.native import resolver + + return resolver.resolve(download=download) + except Exception as exc: # network, unsupported platform, missing cache + print( + f"WARNING: native engine unavailable ({exc}); falling back to " + f"{LEGACY_BINARY}, which cannot be seeded. Reference runs will not " + "be reproducible and single-run comparisons against them are not " + "controlled (issue #228)." + ) + return LEGACY_BINARY + + def run_fortran_amica( data: np.ndarray, params: Dict, @@ -190,12 +217,12 @@ def run_fortran_amica( """Run the AMICA reference binary and collect results. ``binary_path`` selects which reference binary to run. It defaults to the - bundled macOS x86_64 fixture (``pamica/sample_data/amica15mac``); pass the - cross-platform native-engine binary (see ``--native-engine`` in ``main``) to - run the real Fortran reference on Linux/Windows/Apple-Silicon instead. + native engine (see :func:`default_reference_binary`), which is seedable and + therefore reproducible; pass ``LEGACY_BINARY`` for the bundled macOS x86_64 + fixture instead. """ if binary_path is None: - binary_path = Path("pamica/sample_data/amica15mac") + binary_path = default_reference_binary() # Resolve to an absolute path now, before the run chdirs into fortran_dir. binary_path = Path(binary_path).resolve() if not binary_path.exists(): @@ -260,6 +287,17 @@ def run_fortran_amica( print(f"Stdout: {result.stdout}") return None + # A binary predating sccn/amica PR #54 has no `seed` case and its parser + # has no `case default`, so it ignores the keyword in silence and + # re-randomizes instead. Detect that rather than reporting a comparison + # against a random draw as if it were controlled (issue #228). + if "seed =" not in result.stdout: + print( + f"WARNING: {binary_path} did not acknowledge the seed, so it " + "predates the seedable build. Its initialization is random per " + "run and this comparison is not controlled (issue #228)." + ) + except subprocess.TimeoutExpired: os.chdir(original_dir) print("Fortran AMICA timed out")