Seed the Fortran reference and default to the native engine

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 ''}"

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

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

pamica/tests/test_fortran_param_forwarding.py

@@ -105,3 +105,54 @@ 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"
+
+
+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()

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(
@@ -160,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,
@@ -170,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():
@@ -200,11 +247,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
@@ -232,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")