---
title: "Reproducibility, parity and the materials datasheet"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Reproducibility, parity and the materials datasheet}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = FALSE,
  comment = "",
  R.options = list(
    cli.num_colors = 1,
    cli.hyperlink = FALSE,
    crayon.enabled = FALSE,
    width = 80
  )
)
# Console colour carries no meaning on a rendered page. pkgdown turns it on for
# its own build, and the escape sequences then reach the reader as literal text,
# so colour is switched off here for a plain vignette render and a site build
# alike. The fixed width keeps printed output inside the documentation column.
```

A stimulus set is usually described rather than shared. The methods section says
the words were matched on length and frequency, and the set itself, the criteria
that produced it and the code that applied them are not recoverable from the
paper. `lexsync` treats the set as an artefact in its own right: the selection
and the trial order are deterministic, the R and Python engines produce
the same artefacts from the same input, and each run emits a provenance record
with the checksums needed to verify it.

This vignette explains the guarantee and its limits, the mechanisms that achieve
it, the continuous-integration test that keeps it honest, and the two records a
run leaves behind.

```{r setup}
library(lexsync)
schema <- yaml::read_yaml(
  system.file("extdata", "schema.yaml", package = "lexsync")
)
lex <- load_lexicon(
  system.file("extdata", "en_example.csv", package = "lexsync"),
  schema, language = "english"
)
pool <- build_pool(lex, list(length = c(3, 7), frequency = c(3.8, 7.0)))

design <- list(
  name = "vignette_repro", language = "english", n_per_condition = 15,
  conditions = list(
    list(name = "high", define_by = list(frequency = c(5.2, 7.0))),
    list(name = "low",  define_by = list(frequency = c(3.8, 4.4)))
  ),
  match_on = list("length", "n_density", "old20"),
  counterbalance = list(lists = 1)
)
```

## What is guaranteed

The claim is worth stating exactly. Given the same lexicon and the same design,
the R and Python engines produce byte-identical stimuli under the
`standardised_euclidean` and `joint` matching methods, under continuous selection,
and under both pseudoword generators. The claim covers which items are chosen, how
they are paired and the order they are presented in, so the generated experiments
are byte-identical in full. The only exceptions are the `mahalanobis` and
`optimal` matching methods.

'Byte-identical' is meant literally, and stating it that way turned out to matter. The
parity test originally read both engines' stimuli files back with a CSV parser and
compared the values, under which `1` and `1.0` are the same number. On that test 13 of
the shipped designs passed while differing byte for byte. Three of the differences were
serialisation: a whole number written `1` by one writer and `1.0` by the other, a boolean
written `FALSE` and `False`, and a small value written `9e-4` and `0.0009`. One was not:
two reported means differed in the last decimal published, because one engine's
summation was pairwise and the other's was not. Every sum in the package now goes through
a single compensated-summation algorithm written out in both engines, whose agreement
follows from IEEE-754 requiring addition and subtraction to be correctly rounded, and the
mean, variance and standard deviation are built on it. The median sorts and takes the
middle value, or the plain-double average of the two middle values when the count is
even, which keeps it clear of the long-double accumulator that `stats::median()` reaches
by way of `mean()`. The artefacts are now compared as bytes.

Two files sit outside the claim by design, because each records the engine that produced
it: the materials datasheet names its engine and that engine's package versions, and the
run log adds wall-clock timestamps. Their contents are compared field by field instead.

Within one engine, selection is a pure function of its input.

```{r determinism}
a <- match_stimuli(pool, design, schema)
b <- match_stimuli(pool, design, schema)
identical(a, b)
```

That is not a small claim to make about an R function. Repeated calls agree
because nothing in the selection path consults a random number generator, and
because every operation that could depend on the ambient environment has been
pinned. The interesting question is why the *Python* engine also agrees, since
the two implementations share no code.

## How byte-identical output is achieved

Cross-engine identity is not an accident of two careful translations. Each
mechanism below closes a specific way the engines could otherwise drift.

The first is that there is no random number generator anywhere. The obvious way
to choose 15 items from 700 is to sample them, and `lexsync` never does. The
anchor condition is chosen by an even spread across its sorted subpool, taking
positions along the whole range, and every later item is determined by a
distance. A seeded RNG would not help here, because R's Mersenne-Twister stream
and NumPy's are different streams: the same seed gives different draws, so any
sampling step would break parity by construction. Determinism through absence is
the only option that survives two languages, and it now extends to trial order,
whose shuffle is described below.

```{r even-spread}
anchor_pool <- pool[pool$frequency >= 5.2 & pool$frequency <= 7.0, ]
anchor_pool <- anchor_pool[
  order(anchor_pool$frequency, anchor_pool$word, method = "radix"),
]
idx <- unique(round(seq(1, nrow(anchor_pool), length.out = 15)))
idx
```

The spread is `seq()` over the row positions, rounded to integers. Both engines
compute the same positions from the same row count, and the rounding rule is
written out in both.

Sorting is by byte order everywhere. R's default string sort is locale-collated,
and a French locale, a C locale and a Windows locale can order the same words
differently, while any of those orders would differ from Python's, which sorts by
code point. Every sort on the selection path therefore passes `method = "radix"`,
which sorts by bytes.

```{r radix}
words <- c("zebra", "Apple", "apple", "Zebra")
sort(words, method = "radix")
```

Uppercase sorts before lowercase because that is the byte order, and this is what
Python's `sorted()` produces. A locale-collated sort would typically interleave
them. Case folding is pinned for the same reason: `load_lexicon()` lower-cases
through ICU at the root locale rather than through `tolower()`, whose behaviour
depends on the C library and the locale, and which differs from Python's on
Greek final sigma and the dotted Turkish capital I.

Frequencies of bigrams and of subsyllabic constituents are accumulated as
integers, never as running floating-point sums, so the counts cannot depend on
summation order.

Rounding comes before the tie-break, which is the subtle one of these mechanisms.
Two candidate words can sit at exactly equal distance from an anchor, and the
last bits of a square root can differ between two linear-algebra implementations.
If the engines compared raw doubles, an inconsequential difference in the
fifteenth decimal would make one engine call a tie and the other not, and they
would pick different words. Distances are therefore rounded to nine decimal
places before ranking, which puts both engines on the same side of every
comparison that matters.

Rounding creates ties rather than removing them, so the ties in turn need a rule
that both engines can apply identically, and that rule is a total order.
Candidates are ordered by distance, then by the word's bytes, then by its integer
`id`. That is the comparison the matcher makes, and it is reproduced here on
three imaginary candidates.

```{r tie-break}
distance <- c(0.5, 0.5, 0.2)
word     <- c("beta", "alpha", "gamma")
id       <- c(1L, 2L, 3L)
word[order(distance, word, id, method = "radix")]
```

`gamma` wins on distance. `alpha` and `beta` are tied on distance to the ninth
decimal, and the word's bytes settle it, identically in both engines. Since `word`
is unique within a lexicon the order is total, so there is always exactly one
winner and never a coin to flip. The schema notes that this order is fixed in the
code and deliberately not configurable, because it is the mechanism of parity.

The even-spread rounding is shared as well. `unique(round(seq(...)))` is used
identically in both engines, including the `unique()`, which matters when the
pool is smaller than the requested n and the rounding produces a repeated index.

### Trial order: a seeded shuffle without a generator

Trial order used to be the honest exception. Each engine drew the order from its
own seeded generator, R's Mersenne-Twister against NumPy's PCG64, and the same
seed could not give the same permutation, so an order was reproducible within an
engine and different between the engines. That exception is retired. Order is now
decided by a keyed hash: each row of a list is ranked by the SHA-256 digest of
`seed|replicate|list|set|condition`, a tuple that identifies the trial uniquely
under either counterbalancing recipe. Distinct inputs to SHA-256 behave as
independent uniform draws, so sorting rows by their digests realises a seeded
random permutation, and because a digest is a pure function of its input, the
permutation is a pure function of the design: the same bytes from both engines on
any platform, no generator state to save or restore, and a different order for
every seed. The implementation is `.shuffle_deterministic()` in
`R/counterbalancing.R`, mirrored in the Python engine's `counterbalancing.py`.

```{r trial-order}
stim <- counterbalance(match_stimuli(pool, design, schema), design, schema)
head(stim[, c("trial", "condition", "word")], 4)
same <- counterbalance(match_stimuli(pool, design, schema), design, schema)
identical(stim$trial, same$trial)
```

The order is still randomised in the sense a study needs: an item's properties
bear no systematic relation to its position, and a different seed gives a
different order. What changed is that the permutation is engine-invariant, so it
now sits inside the cross-engine contract, and the guarantee reaches the
generated experiments whole: for every shipped design, the PsychoPy script, the
OpenSesame file, the jsPsych page and both loop-table CSVs are byte-identical
across the engines.

One detail matters to anyone embedding `lexsync` in a longer script:
counterbalancing no longer touches the random number generator at all, since
nothing in the package consults one, so calling it cannot disturb the stream your
own code is drawing from.

## The Unicode input contract

Word matching, deduplication and the `length` dimension operate on Unicode code
points exactly as supplied. Beyond the pinned lower-casing described above, the
orthographic form is never rewritten, and in particular no canonical
normalisation is applied, neither NFC nor NFD. A precomposed character such as
`é` and its decomposed counterpart, the base letter followed by a combining
accent, are therefore distinct items that differ in length, and a lexicon
mixing the two conventions carries apparent duplicates that never match one
another. Supply NFC-normalised text, which is the form most corpora already
distribute.

The syllable estimator and the pseudoword generators carry a narrower contract
still. The vowel class behind `n_syllables` covers the Latin-1 vowels, and the
pseudoword generators build their candidates from lower-case a to z, so both
are orthographic approximations defined for Latin orthographies and unavailable
for other scripts. The shipped Chinese design matches on `length`, `n_density`
and `old20` for exactly this reason.

## The documented exceptions

`mahalanobis` and `optimal` are equivalent to their byte-identical siblings in
what they compute and not in what they select, and the reason is worth
understanding.

`mahalanobis` needs the inverse of the pool's correlation matrix. A matrix
inverse is computed by LAPACK, and R's LAPACK is not necessarily the one NumPy is
built against. The two may differ in blocking, in vectorisation and in the order
they accumulate sums. The results agree to well within any tolerance you would
care about and they need not agree in their last bits, and the rounding that saves
the Euclidean path cannot save this one: the distance is now a function of a
matrix whose entries already differ below the rounding threshold, so the rounding
can land on either side.

`optimal` solves a linear assignment problem. Where several assignments have the
same total cost, which one a solver returns is an implementation detail of the
solver, and R's `clue::solve_LSAP()` and SciPy's implementation resolve those ties
differently. Both return an optimal assignment. They need not return the
same one.

Neither is a bug that a determined engineer could fix. Guaranteeing them would
mean shipping a bit-exact matrix inverse and a tie-canonical assignment solver in
both languages, which is a large amount of numerical code to maintain for a modest
gain. The chosen trade-off is to keep the byte-identical methods as the defaults,
document the exception, and record it per run.

```{r cross-engine}
maha <- design
maha$matching <- list(method = "mahalanobis")
maha_stim <- match_stimuli(pool, maha, schema)
maha_ds <- build_datasheet(
  maha, schema, NULL, maha_stim,
  system.file("extdata", "en_example.csv", package = "lexsync"),
  list(stimuli = NA_character_), schema$seed
)
maha_ds$selection$cross_engine
```

Every datasheet carries this field, so a reader of the record does not have to
know the rule. It reads `byte-identical` for the default methods and
`approximate (platform linear algebra)` for these two, and the generated methods
paragraph adds a sentence stating the caveat in prose when it applies. The
constraint also shapes the roadmap: a covariance-aware distance might well be
the better default, and it stays optional precisely because it cannot carry the
guarantee.

One further place is worth naming, though it is not an exception. Almost
everything on the path to a compared artefact is built from operations IEEE-754
either mandates correctly rounded or makes exact, so the engines agree by
construction rather than by measurement. The Student t quantile and distribution
behind the confidence interval and the TOST *p*-value are where that argument
does not reach: `stats::qt` and `stats::pt` on one side, SciPy's on the other, two
independent implementations that were measured to disagree by up to about 1e-12
over the range these reports use. The published values are rounded to three and
four decimal places, which leaves 5e-4 and 5e-5 of headroom, so the comparisons
files still come out byte-identical. The guarantee there rests on a margin of some
seven orders of magnitude rather than on the construction, and saying so is what
stops anyone widening the reported precision without checking the margin first.

## The parity test

A guarantee that is not tested is a comment. The parity gate in continuous
integration runs the R pipeline over the shipped designs, regenerates each with
the Python engine and compares the two, across designs covering both engines'
matching methods, continuous selection, both pseudoword generators, resampling,
priming and self-paced reading, three languages including a logographic script,
and three reproductions of published designs.

The comparison is on identity columns: the word or target, the condition and the
`set` index. Including `set` is what gives the test teeth, since without it the
engines could select the same items and pair them differently and still pass.
Having matched on those, the test then checks every other shared column, so the
`trial` position and the dimensions each engine computes at run time
(`n_density`, `old20`, `n_syllables`, `bigram_freq`) must also agree.

Three details keep the gate from passing vacuously. The test skips when the R
reference output is absent, since the Python package can be installed on its
own, so CI sets `LEXSYNC_REQUIRE_PARITY=1` to turn those skips into failures: a
job meant to run both engines cannot pass by quietly skipping. The R references
are regenerated and then checked with `git diff --exit-code`, so an R-engine
regression fails the build. Without that step the comparison would run against a
stale snapshot that the regression matches. And the R suite is run once from the
source tree with a check that nothing skipped, because the repository-coupled
tests silently skip when run from an installed copy.

## The materials datasheet

Each run emits a datasheet in two forms, a JSON record for machines and a
Markdown rendering for people. It exists because the reproducibility literature
in linguistics keeps finding the same gap: materials are described but rarely
deposited in a form that permits reuse or verification (Bochynska et al., 2023),
and the degrees of freedom in preparing them are large enough to matter
(Roettger, 2019). A datasheet is the FAIR-style record of one selection
(Wilkinson et al., 2016).

```{r datasheet}
out <- file.path(tempdir(), "lexsync_repro")
dir.create(out, showWarnings = FALSE)
report <- match_report(
  stim, c("frequency", "length", "n_density", "old20"), schema
)
stim_path <- file.path(out, "stimuli.csv")
write.csv(stim, stim_path, row.names = FALSE)

ds <- build_datasheet(
  design, schema, report, stim,
  source_path = system.file("extdata", "en_example.csv", package = "lexsync"),
  artifacts = list(stimuli = stim_path),
  seed = schema$seed,
  candidate_pool = lapply(design$conditions, function(cnd)
    list(
      condition = cnd$name,
      n_candidates = nrow(build_pool(pool, cnd$define_by))
    ))
)
names(ds)
```

The record answers the questions a sceptical reader would ask. Where the words
came from, and whether this is the same file, is settled by `materials_source`,
which carries the path and its SHA-256. What was applied is recorded in
`selection`, covering the method, the matched dimensions and the tolerance
windows as resolved, which is the design's overrides merged onto
the schema's defaults rather than the defaults a design may have replaced. What
the selection achieved is in `realised_control`, as the effect sizes, intervals,
variance ratios and TOST verdicts. What it could have chosen instead is in
`candidate_pool`, which records how many items satisfied each condition before
matching.

```{r checksums}
ds$materials_source$sha256
str(ds$selection$tolerance_k)
do.call(rbind, lapply(ds$selection$candidate_pool, as.data.frame))
```

The candidate-pool counts do a job that is easy to miss. Here the high condition
was drawn from a few dozen eligible words and the low condition from several
hundred, which are different epistemic situations: choosing 15 items from a large
candidate set invites the concern that a selection could be steered towards almost
any result, whereas a condition with barely more candidates than items had little
freedom to be steered. Recording the counts, alongside the fact that the selection
is deterministic and blind to any outcome measure, lets a reader judge that rather
than take it on trust.

Note that the resolved windows include a `frequency` entry even though this design
matches on the three orthographic dimensions only. The record reports the windows
as resolved from the schema and the design, and a window is applied only to a
dimension named in `match_on`.

```{r versions}
str(ds$reproducibility)
```

The reproducibility block records the seed and the versions of the engine and its
dependencies, since a bug fix in a dependency can change a selection and the
record should say which versions produced this one.

### The methods paragraph and the pre-registration template

The datasheet renders itself into prose, which is a small feature that removes a
common source of error: a methods section transcribed by hand from numbers on a
screen.

```{r methods-paragraph}
cat(strwrap(methods_paragraph(ds), width = 76), sep = "\n")
```

The paragraph is generated from the record, so it cannot disagree with the
stimuli. It states the candidate-pool size, and it adds the cross-engine caveat
when the design's method warrants it. On the realised control it takes one of two
branches. Where every controlled dimension passes its equivalence test, it names
the largest standardised difference on any matched dimension with its interval,
rather than a selective best case, and says that this sits within the equivalence
bound. Where one does not, it declines to summarise and sends the reader to the
realised-control table instead, which is the branch this run takes.

That is worth seeing, because the paragraph is a draft to adapt rather than one to
paste unread. What fails here is not the matching but the width of the interval,
and the interval is wide because this demonstration matches only 15 items per
condition, which is precisely the dependence on set size that the matching
vignette warns against reading past. A real study would either match enough items
to narrow the interval or say plainly that the control is established at the point
estimate and not at the bound.

`write_datasheet()` writes both forms, and the Markdown rendering carries a
pre-registration skeleton with its materials section already filled in.

```{r write-datasheet}
paths <- write_datasheet(
  ds, file.path(out, "datasheet.json"), file.path(out, "datasheet.md")
)
md <- readLines(file.path(out, "datasheet.md"))
grep("^#", md, value = TRUE)
```

The sections left blank are the ones no tool can fill: the hypotheses, the
sampling plan and the inference criteria. The analysis plan is the exception: it
arrives with a suggested model, and the suggestion takes a position.

```{r analysis}
ds$analysis$suggested_model
```

The model crosses random effects for subjects and items because items are a sample
of the language rather than the population of interest, and an analysis treating
them as fixed generalises only to the words that were used (Clark, 1973; Yarkoni,
2022). The accompanying note recommends beginning with the maximal structure and
reducing it if it does not converge. It is a starting point for a design of this
shape.

## The run log

The datasheet records what a run produced. The run log records what it did.

```{r run-log}
log <- new_run_log(
  "vignette_repro",
  meta = list(seed = schema$seed, language = "english")
)
log <- log_step(
  log,
  sprintf("lexicon loaded: %d words", nrow(lex)),
  list(words = nrow(lex))
)
log <- log_step(log, sprintf("pool after filters: %d words", nrow(pool)))
log <- log_artefact(log, stim_path, rows = nrow(stim))
log_md <- write_run_log(log, file.path(out, "run_log.md"),
                        file.path(out, "run_log.jsonl"))
cat(readLines(log_md), sep = "\n")
```

The pipeline builds this log as it goes, recording each stage with its parameters,
the equivalence verdict on every dimension, any balance warning, and an entry per
file written with its row count and MD5. Both a Markdown log and a JSON Lines log
are emitted, the latter so that a run can be checked by a script. The MD5 here is a
content fingerprint for spotting drift between runs, whereas the datasheet's
SHA-256 is the stronger one that belongs in a permanent record.

## Reproducing a run

The pipeline does all of the above in one call, and `run_pipeline()` on a design
configuration is the reproducible unit: it writes the stimuli, the reports, the
three experiments, the datasheet and the log. Sharing the design file and the
corpus is then enough for someone else to regenerate the materials and verify the
checksums, in either language.

Two habits make that work. Pin the corpus, since a selection is only reproducible
against a fixed lexicon, which is why the bundled corpora are a checksummed
snapshot and why the `wordfreq` connector draws on a package frozen in 2024, whose
data is a stable snapshot of usage through roughly 2021. And where cross-engine
reproducibility matters, prefer the default matching methods, since they are the
ones that carry the guarantee.

## References

Bochynska, A., Keeble, L., Halfacre, C., Casillas, J. V., Champagne, I.-A., Chen,
K., Röthlisberger, M., Buchanan, E. M., & Roettger, T. B. (2023). Reproducible
research practices and transparency across linguistics. *Glossa
Psycholinguistics*, *2*(1). https://doi.org/10.5070/G6011239

Clark, H. H. (1973). The language-as-fixed-effect fallacy: A critique of language
statistics in psychological research. *Journal of Verbal Learning and Verbal
Behavior*, *12*(4), 335–359. https://doi.org/10.1016/S0022-5371(73)80014-3

Roettger, T. B. (2019). Researcher degrees of freedom in phonetic research.
*Laboratory Phonology*, *10*(1), Article 1.
https://doi.org/10.5334/labphon.147

Wilkinson, M. D., Dumontier, M., Aalbersberg, IJ. J., Appleton, G., Axton, M.,
Baak, A., Blomberg, N., Boiten, J.-W., da Silva Santos, L. B., Bourne, P. E.,
Bouwman, J., Brookes, A. J., Clark, T., Crosas, M., Dillo, I., Dumon, O.,
Edmunds, S., Evelo, C. T., Finkers, R., … Mons, B. (2016). The FAIR Guiding
Principles for scientific data management and stewardship. *Scientific Data*,
*3*, Article 160018. https://doi.org/10.1038/sdata.2016.18

Yarkoni, T. (2022). The generalizability crisis. *Behavioral and Brain Sciences*,
*45*, Article e1. https://doi.org/10.1017/S0140525X20001685