Add a script that reproduces the benchmarks on any machine

CHANGELOG.md

@@ -9,6 +9,12 @@ the README for what that means before 1.0.
 
 ### Added
 
+- `bench/ramtools-bench.py` reproduces the benchmarks on any Linux or macOS
+  machine: it downloads chr20 of NA12878 from EBI and checks it against known
+  md5 sums, builds two revisions the same way, measures conversion, size, the
+  round trip, dumps, regions and statistics in alternating order, checks every
+  output against samtools, and writes a report. See `bench/README.md`
+  ([#116]).
 - `-fqzcomp` on both converters compresses lossless quality scores with
   fqzcomp, the codec of CRAM 3.1, in blocks of 10,000 records. On 2 million
   records of chromosome 11 the file is 8.1% smaller, conversion takes 20 to 30%
@@ -264,3 +270,4 @@ it stay readable by every later release.
 [#113]: https://github.com/compiler-research/ramtools/pull/113
 [#114]: https://github.com/compiler-research/ramtools/pull/114
 [#115]: https://github.com/compiler-research/ramtools/pull/115
+[#116]: https://github.com/compiler-research/ramtools/pull/116

CONTRIBUTING.md

@@ -87,7 +87,10 @@ Every pull request runs four checks:
   the summary.
 - Put the evidence in the description: what you ran, what it printed, and what
   you compared it against. Reviewers check claims against samtools, so give
-  them the command that shows the agreement.
+  them the command that shows the agreement. For a change to speed, memory or
+  size, attach the report of `bench/ramtools-bench.py run` (see
+  `bench/README.md`) on 2 million records, and on all of chr20 when the file
+  format changes.
 - Add a line under `[Unreleased]` in `CHANGELOG.md`, in Added, Changed, Removed
   or Fixed, with a link to the PR. Say what a user notices and, for speed or
   size, by how much. Changes users do not see, such as tests, CI or a

README.md

@@ -212,6 +212,11 @@ The two counts agree, the `diff`s print nothing, and for a SAM input
 
 ## Benchmarks
 
+The numbers in this section predate the current file format.
+[bench/README.md](bench/README.md) measures the current code on your
+own machine: it downloads chr20 of NA12878, builds two revisions, measures
+them against each other and against samtools, and writes a report.
+
 HG00154 from the 1000 Genomes Project: 196,040,370 records, 72.06 GB as SAM,
 coordinate sorted, GRCh37 reference names (`1`, not `chr1`). Run on an Intel
 i3-7020U laptop (2 cores, 4 threads, 12 GB of memory) with input and output on

bench/README.md

@@ -0,0 +1,120 @@
+# Benchmarks on real data
+
+`bench/ramtools-bench.py` reproduces the sizes and timings quoted in
+`CHANGELOG.md` on any Linux or macOS machine. It downloads the data, builds two
+revisions of ramtools the same way, measures them in alternating order, checks
+every output against samtools, and writes a report.
+
+## What you need
+
+- Linux or macOS, Python 3.8 or newer, git, CMake 3.16 or newer and a C++17
+  compiler.
+- ROOT 6.38 or newer, set up in the shell (`source <ROOT>/bin/thisroot.sh`).
+- htslib with its development files, as for building ramtools.
+- samtools 1.16 or newer, built with libcurl, which reads the data from EBI.
+- About 15 GB of disk for the data, and 3 GB more while all of chr20 is
+  measured.
+
+`bench/ramtools-bench.py check` lists what is there and what is missing.
+
+## Running it
+
+```bash
+git clone https://github.com/compiler-research/ramtools.git
+cd ramtools
+source <ROOT>/bin/thisroot.sh
+bench/ramtools-bench.py check
+bench/ramtools-bench.py fetch                                  # once: 20 to 60 minutes
+bench/ramtools-bench.py run --base v0.1.0 --head develop       # 2 million records
+bench/ramtools-bench.py run --base develop~1 --head develop --input chr20
+```
+
+`--base` and `--head` take anything git does: a branch, a tag, a commit. The
+default head, `WORKTREE`, is the checkout with its uncommitted changes. To
+measure a pull request, fetch it first:
+`git fetch origin pull/115/head:pr-115`, then `--head pr-115`.
+
+Everything goes into `--work`, by default `../ramtools-bench-work` next to the
+clone, or `$RAMTOOLS_BENCH_WORK`: the data, a build per revision (built once
+and kept), and the results.
+
+On a machine with several processor sockets, `--pin 0` runs every measured
+command on NUMA node 0 with `numactl`, so that threads and memory stay on one
+socket. On Linux, the `performance` CPU frequency governor gives steadier
+timings; `check` says which governor is set.
+
+## The data
+
+| | |
+|---|---|
+| Sample | NA12878, 1000 Genomes 30x high-coverage set (NYGC), GRCh38, NovaSeq, qualities binned to four values |
+| Source | `https://ftp.sra.ebi.ac.uk/vol1/run/ERR323/ERR3239334/NA12878.final.cram`, chr20 only |
+| Reference | chr20 of `https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa`; md5 of the sequence `b18e6c531b0bd70e949a7fc20859cb01`, the M5 of chr20 in the CRAM |
+| Records | 17,705,654; md5 of `samtools view` of the BAM, without the header, `89386f84300ce2a3b5100ec5760a197f` |
+
+`fetch` checks the reference and the records against these and stops when
+they differ; `--accept-different-data` goes on, and the report then shows
+another md5 for the SAM. It writes, in `--work`:
+
+- `GRCh38.chr20.fa`: `samtools faidx <reference> chr20`.
+- `NA12878.chr20.bam`: `samtools view --no-PG -b -T GRCh38.chr20.fa -X <CRAM> <CRAI> chr20`, indexed.
+- `NA12878.chr20.sam`: `samtools view --no-PG -h` of the BAM, about 9 GB.
+- `NA12878.chr20.cram`: CRAM 3.1 with samtools' default profile, indexed, for
+  the comparison with samtools.
+- The first 2,000,000 and 50,000 records of the SAM, each also as BAM and CRAM.
+
+`--no-PG` keeps the commands out of the header, so every machine gets the same
+files: the report gives the md5 of the SAM it measured.
+
+## What is measured
+
+`run --input` picks the records: `small` (50,000, a minute or two, to try the
+script out), `subset` (2 million, the default, for comparing changes) or
+`chr20` (all 17.7 million, once per change). The steps:
+
+| Step | Command | small | subset | chr20 |
+|---|---|---:|---:|---:|
+| convert SAM, 1 thread | `samtoramntuple in.sam out.ram` | 3 | 5 | 0 |
+| convert SAM, threads | `samtoramntuple in.sam out.ram -threads 4` | 3 | 5 | 2 |
+| convert BAM | `bamtoramntuple in.bam out.ram` | 3 | 5 | 1 |
+| whole file | `ramdump -h` | 3 | 5 | 3 |
+| large region | `ramdump file region`, about half the records (10 Mb on chr20) | 3 | 5 | 3 |
+| three small regions | `ramdump file region`, 500 to 2,000 bp | 5 | 5 | 3 |
+| count | `ramdump -c` of the large region | 5 | 5 | 3 |
+| statistics | `ramstats flagstat`, `ramstats idxstats` | 5 | 5 | 3 |
+
+The numbers are the timed runs per revision; `--repeat` and `--repeat-convert`
+change them. The two revisions take turns, base and head, then head and base,
+so a machine that slows down part of the way slows both alike. Reads have one
+untimed run first, so the file is in the page cache for every timed run.
+
+For each run the script records the wall time, the CPU time (user plus
+system, of all the process's threads), the peak memory, and the load average
+before it. Peak memory is the peak physical footprint on macOS and the peak
+RSS on Linux; the two are not comparable with each other.
+
+samtools runs the same reads on the BAM and on the CRAM 3.1, one thread each,
+as the yardstick (`--no-samtools` leaves it out).
+
+After the timings, every result is checked: `ramdump -h` of the file converted
+from SAM gives the SAM byte for byte; `ramdump` of the file converted from BAM
+gives what `samtools view` prints; the count, the small regions, `flagstat`
+and `idxstats` equal samtools' on the BAM. `run` exits with 1 when a check
+fails.
+
+## Reading the results
+
+`run` prints the report and saves it as Markdown and JSON in `WORK/results`,
+named after the time, the machine, the input and the two builds. The JSON has
+every run; `bench/ramtools-bench.py report <file.json>` prints the report
+again.
+
+Each time is the median of the runs, with the fastest and the slowest in
+brackets. Compare CPU times when the machine is busy: wall times then vary
+with the load, CPU times much less. When the brackets of base and head
+overlap, the difference is not shown by this run; run more repeats, or run
+again when the machine is quiet. A size, a correctness check and an output are
+exact.
+
+To compare two machines, run `fetch` on both and check that the reports give
+the same md5 for the SAM.