Reproducibility

A number in a paper should be something you can produce again: next month, on another machine, or by a referee who has your script. That needs four things, and Mera supports each one.

What this page assumes

That you can start Julia and install a package. If package mode (the ] prompt) is new to you, read Julia for Simulation Analysis first: it introduces environments in a few lines, and this page is the full treatment of the same idea.

Set this up before you start analysing, not after. Reconstructing which versions produced a figure six months ago is difficult; recording it as you go costs nothing.

FAIR in 5 lines

FAIR stands for Findable, Accessible, Interoperable and Reusable: results that others can find, open with their own tools, and repeat. Mera does most of the work: every file it saves gets a sidecar, a small <file>.meta.json next to it that records which snapshot, code and package versions made the file.

Start each analysis in its own folder: meraproject("galaxy_paper") makes it, then start Julia there with julia --project=. (the first time: ] instantiate). The folder keeps its own list of package versions, so the sidecars can record exactly what ran. Then:

exportdata(m, "fig3.h5")                            # 1. saves write their sidecar by themselves
save("fig3.png", fig); writesidecar("fig3.png", m)  # 2. a figure: one more line
readsidecar("fig3.png")                             # 3. where a file came from
sidecaroverview(".")                                # 4. all files of the analysis
bundle("."; simulation=info)                        # 5. before publishing: one record, and the BibTeX

Nothing else has to be switched on. If you skip the first step, Mera reminds you once. The rest of this page explains each part.

what you needhow
the same package versionsa Julia project, described on this page
a record of what produced a numberprovenance on an object, a sidecar next to each saved file, bundle for a folder: Sharing Results & Sidecars
the formula behind a derived quantityHow Quantities Are Computed
confidence the package itself is correcttest simulations with known answers

The first is the one people skip, so it comes first here.

Use a Julia project

By default, packages you install go into a global environment shared by everything on your machine. Installing an unrelated package for a different task can then change the versions used by your analysis, without you touching it. A project avoids that: it is a directory with its own package list, so each analysis has its own fixed set of versions.

Create one for your analysis:

julia> ]                        # press ] to enter package mode
(@v1.12) pkg> activate .        # this directory is now the project
(my_analysis) pkg> add Mera CairoMakie

The prompt changes from @v1.12 to my_analysis, which tells you the project is active. Two files appear in the directory:

my_analysis/
├── Project.toml     # what you asked for:  Mera, CairoMakie
├── Manifest.toml    # what you actually got: exact versions of those AND everything they depend on
└── analysis.jl

The difference between them matters:

  • Project.toml lists your direct dependencies, with loose version bounds. It says "some Mera 1.x".
  • Manifest.toml records the exact resolved version of every package in the tree, direct and indirect, hundreds of them. It says "Mera 1.8.0, JLD2 0.6.2, and 300 more, exactly these".

Manifest.toml is what makes a result reproducible. Keep it next to your script and commit it to version control.

Mera can also start the folder for you: meraproject("my_analysis") writes a Project.toml with the Mera you run now, a first analysis.jl, an analysis.toml for the title, authors and license, and a short README. ] instantiate in that folder then installs Mera there and writes the Manifest.toml.

Start Julia in the project

The project is only used if you ask for it:

julia --project=.                # from inside the project directory
julia --project=/path/to/my_analysis

Forgetting this is the most common mistake. Without it you are back in the global environment and may get different versions than the ones you recorded. To check inside a session:

using Pkg
Pkg.status()          # shows the active project path at the top

In VS Code, select the project folder as the Julia environment (bottom status bar). In a Jupyter notebook, activate it in the first cell:

using Pkg; Pkg.activate(@__DIR__)

Reproduce it later, or elsewhere

Given the directory with both TOML files, anyone rebuilds the identical stack:

git clone https://example.org/my_analysis && cd my_analysis
julia --project=. -e 'using Pkg; Pkg.instantiate()'
julia --project=. analysis.jl

instantiate reads Manifest.toml and downloads exactly those versions. It does not resolve anything, so the result is the same package set you had, not merely a compatible one.

Two limits worth knowing:

  • A Manifest.toml is tied to a Julia minor version. One written on 1.11 may not instantiate on 1.10. Record the Julia version too, versioninfo() prints it, and say it in your methods section.
  • Packages with binary dependencies resolve per platform. The Julia packages match; the underlying system libraries may differ between macOS and Linux.

A complete example

# analysis.jl, run with:  julia --project=. analysis.jl
using Mera

info = getinfo(300, "/data/sim/mw_L10")
gas  = gethydro(info, lmax=10)

m = msum(gas, :Msol)
println("total gas mass: ", m, " Msol")

# record what produced that number, next to the number itself
p = provenance(gas)
println(p)

provenance returns the Mera version and build, the simulation code and output, the units in use and a timestamp. Print it into your log, or save it beside the figure. Together with the committed Manifest.toml it answers "what produced this?" completely. See Provenance for the full record and how to store it.

Maps and camera views

provenance describes the snapshot and the build, not the view. A projection map carries its view itself: m.los, m.up, m.center, m.pixsize and m.extent, and savemap keeps them together with the provenance. One thing is not recorded: direction=:faceon, :edgeon and axis=:angmom compute the axis from the data that gets drawn (what you loaded, the mask, the fov sphere), so a different load or selection gives a different camera. For a result you want to repeat, store the vector (angmom, or m.los and m.up) and pass it back in. The off-axis cheatsheet shows the call.

Sidecar files

Every file Mera writes for you, a Mera file from savedata, a map from savemap, a movie from savemovie, a clump catalogue from save_clumps or an export from exportdata, gets a readable record next to it. A figure is saved by the plotting package, so it needs one more line, writesidecar:

exportdata(m, "fig3.h5"; length_unit=:kpc, note="Figure 3")
# writes fig3.h5 and fig3.h5.meta.json

save("fig3.png", fig)                  # Makie
writesidecar("fig3.png", m)            # writes fig3.png.meta.json

The sidecar holds the provenance of the snapshot, the units and geometry of the result, the SHA-256 of the file itself and of the run's namelist.txt, makefile.txt, patches.txt and compilation.txt, the environment that ran (Julia and Mera version, operating system, CPU model, thread counts), and the analysis around it: the SHA-256 of your Project.toml and Manifest.toml, the versions of your project's packages, and the script that ran, by name and SHA-256 (include_script=true also stores its text). Publish it with the data, for example in the same Zenodo record. It is JSON, so a reader needs neither Julia nor Mera.

Two things it does not record. The host name is left out on purpose, and anonymize=true also shortens the simulation path to its last two parts. And a function you wrote yourself, such as the reducer of a timeseries, appears only by name: the script or notebook that ran is the record of what it did, so publish it next to the results.

To read one back, readsidecar("fig3.h5") prints it and says whether fig3.h5 still matches the SHA-256 it was written with; sidecaroverview("results/") lists every sidecar in a folder. sidecar=false turns them off for one call; sidecar(false) for the session.

Why they are on by default. A record helps only if it exists when someone asks for it, often a year later, and a feature that has to be switched on is mostly never switched on. A sidecar costs little: about 5 KB, one short line the first time a session writes one, and it never breaks a save, since a sidecar that cannot be written becomes a warning. Switch them off with sidecar(false) at the top of a script when the number of files matters, for example on a cluster with a file quota, or without touching the script: MERA_SIDECAR=0 julia analysis.jl. Use anonymize=true when the simulation path should not be public.

Two set-ups make a sidecar less useful, and Mera says so the first time it meets them in a session. If Julia runs in the shared environment (the @v1.12 prompt) instead of the analysis's own project, the recorded Manifest.toml changes whenever you add any package, so it no longer says which versions made the result. If Mera runs from a checkout with uncommitted changes, the recorded commit is not the code that ran. Sharing Results & Sidecars walks through all of it on a real simulation.

Threads

The number of threads Julia runs with can change the last digits of some results, never more. Measured at 1, 2, 4 and 8 threads, and checked by the test suite:

resultsame thread count, run againdifferent thread count
reading: gethydro, getparticles, ...identicalidentical
axis-aligned hydro projectionidenticalidentical
clumpfindidenticalidentical
off-axis projection, particle projectionidenticalround-off: below 10⁻¹⁴ relative, 5·10⁻¹² for a dispersion
volumerenderidenticalround-off: 10⁻¹⁷ in the colours

So a result repeats exactly when it is computed again with the same number of threads. The off-axis, particle and volume renderings split the cells into one chunk per thread and add the partial images; with a different thread count the additions happen in a different order, which moves the last bits. Every sidecar records the thread count (julia_threads), so the condition for an exact repeat is written down with the result.

Checking Mera itself

The versions being pinned is only useful if the package is right. Mera publishes fifteen small test simulations, each with a known answer, either a law that follows from its own setup or reference values published by the RAMSES developers:

ENV["MERA_TEST_DATA"] = download_testdata()
ENV["MERA_SMOKE_ONLY"] = "1"      # the tier the public simulations can back
using Pkg; Pkg.test("Mera")

MERA_SMOKE_ONLY matters here. Setting MERA_TEST_DATA alone tells the suite that simulation data is present, and it then also reaches for runs that are not published, so it would fail on something you never had. With the flag you get every test the public simulations can actually support.

What each tier of that suite proves, and which simulation backs which test, is set out in the Testing Framework.

That runs the suite against real RAMSES output on your machine, with your versions. Without the data it still runs and passes, covering every analytic correctness check; see test/README.md.

Checking a formula

Every derived quantity is written out in How Quantities Are Computed, transcribed from the source, so you can confirm what a number means rather than trusting a name. If a definition there disagrees with the one your field uses, you will see it immediately, which is the point.

Periodic boxes

RAMSES does not record its boundary conditions in info_*.txt. When a namelist is present, Mera infers them and says so in the getinfo summary:

boundaries:       periodic in x, y, z

info.boundaries is :periodic, :nonperiodic, :mixed, or :unknown when there is no namelist to read.

:mixed is common, not exotic. RAMSES closes individual faces, so a run can wrap in some directions and not others: RAMSES's own rad_beams.nml closes x and y and leaves z periodic. Ask which:

Mera.periodic_axes(info.namelist_content)   # (x = false, y = false, z = true)

and the summary line says so too:

boundaries:       periodic in z only (x, y closed by &BOUNDARY_PARAMS)

Mera reports this and does nothing else with it. No function wraps coordinates on your behalf, because doing that to a run with outflow or reflecting boundaries would give wrong physics with no warning, and the boundary type cannot be recovered from the snapshot alone. Where wrapping matters, you ask for it:

getvar(gas, :r_sphere_periodic, center=[0., 0., 0.])   # minimum-image radius

Two places to watch on a periodic run, because neither wraps:

  • a subregion touching a face is truncated, not wrapped around
  • center_of_mass of a structure straddling a face lands in the middle of the box, which looks plausible and is wrong

A short checklist

  • [ ] the analysis lives in its own directory with Project.toml and Manifest.toml
  • [ ] both are committed
  • [ ] Julia is started with --project=
  • [ ] the Julia version is recorded
  • [ ] the number of threads is recorded (every sidecar does), for results that must repeat bit for bit
  • [ ] provenance output is saved with each result (the sidecar does this for every file Mera writes)
  • [ ] the scripts or notebooks that made the results are published next to them
  • [ ] every figure file has a sidecar (writesidecar right after the save)
  • [ ] the software is cited at its exact version (citation(gas) gives the BibTeX)
  • [ ] the published folder has its ro-crate-metadata.json (bundle("."))
  • [ ] for camera views, the viewing vectors (m.los, m.up, or the angmom axis) are saved too
  • [ ] the simulation output and its path are named in the script, not typed by hand each run