diff --git a/.gitignore b/.gitignore index 01b10c9b..f349fa3a 100644 --- a/.gitignore +++ b/.gitignore @@ -10,13 +10,23 @@ src/target data-stress *.fasta *.fasta.gz +*.fastq +*.fastq.gz +*.vcf *.zst *.zst.meta *.pb -./**/*.json +**/*.json *.bin *.log *.csv +*.meta +*.pfiv +*.done +*.efh +*.efl +*.bin.idx +*.prsb Betula_exilis--IGA-24-33 benchmark/genomes benchmark/genomes_orig @@ -34,6 +44,7 @@ benchmark/reference_dist benchmark/obikmer_dist benchmark/specific_index_count benchmark/specific_index_presence +benchmark/query_data TNT phyg biblio diff --git a/DevDocMD/implementation/benchmark_query_testing.md b/DevDocMD/implementation/benchmark_query_testing.md index 023bea8a..e234e984 100644 --- a/DevDocMD/implementation/benchmark_query_testing.md +++ b/DevDocMD/implementation/benchmark_query_testing.md @@ -1,84 +1,93 @@ -# Benchmark: query-path testing (discussion) +# Benchmark: query-path testing `benchmark/Makefile` exercises indexing, merge, and phylo distance -reconstruction against simulated bacterial genomes, but has no coverage of -`obikmer query` — the read-matching path — nor of the sparse packed -presence-matrix format (`obikmer pack --sparse`). This note captures the -planned extension. +reconstruction against simulated bacterial genomes. It now also covers +`obikmer query` — the read-matching path — and the sparse packed +presence-matrix format (`obikmer pack --sparse`), previously untested by +this pipeline. ## Motivation -- `query` is untested end-to-end. A regression there would not be caught by - the existing `verify_presence`/`verify_merge_presence` branches, which only - check index *content* against the `.npz` truth, never the query API. +- `query` had no end-to-end coverage. A regression there would not be caught + by `verify_presence`/`verify_merge_presence`, which only check index + *content* against the `.npz` truth, never the query API. - `pack --sparse` produces a presence-matrix format documented (see [siblings.md](../architecture/siblings.md)) as faster for single-row access (query) and slower for column-oriented access (phylo `--metric`). - `global_index_presence/` built by `merge_presence.sh` is currently always - packed dense (packing is a stage inside `merge`, not a separate `pack` - invocation). There is no dense/sparse regression check. + `global_index_presence/` built by `merge_presence.sh` is always packed + dense (packing is a stage inside `merge`, not a separate `pack` + invocation) — there was no dense/sparse regression check. -## Plan +## Query read source -**New read source, independent of `simulated_data/`.** Reusing -`simulated_data///reads_R1.fastq.gz` for queries would bias -the test: those reads were already folded into the index being queried, with -the same sequencing-error draw. Query reads must come from a *second*, -independent `iss generate` run against the same reference genome(s) — new -random error draw, same underlying sequence — landing in a separate tree: -`query_data///reads_R1.fastq.gz`, built by the existing -`simulate_one.sh` (unseeded, so a second invocation naturally draws different -reads). +Query reads are independent of `simulated_data/` (which is folded into the +index being queried): reusing those reads would test against the exact +error draw the index was built from. `query_data///` holds +a *second*, independent `iss generate` run against the same reference +genome, via `simulate_query_one.sh` — unseeded, so a second draw picks up +different sequencing errors than `simulate_one.sh`'s draw for the same +genome. Fixed at 100,000 read pairs per genome (not coverage-proportional +like the 15x used for `simulated_data/`), so wall/RSS numbers stay +comparable across genomes of very different sizes. -Two specimens chosen as query sources (enough to catch a dense/sparse -regression without duplicating the exhaustive per-specimen coverage -`verify_merge_presence` already provides across all `SPECIMENS`): -`Escherichia_coli--K-12_MG1655` (common, well-represented bacterium) and -`Saccharolobus_islandicus--M.16.4` (the only archaeon in `SPECIES` — distant -lineage, different GC content, stresses the query path differently from a -close-relative match). +Two query-source specimens, hardcoded as `QUERY_SPECIMENS` in +`make_deps.py`: `Escherichia_coli--K-12_MG1655` (common, well-represented +bacterium) and `Saccharolobus_islandicus--M.16.4` (the only archaeon in +`SPECIES` — distant lineage, stresses the query path differently from a +close-relative match). Two is enough to catch a dense/sparse regression +without duplicating the exhaustive per-specimen coverage +`verify_merge_presence` already provides across all `SPECIMENS`. -`make_deps.py` needs a `QUERY_SPECIMENS` list (explicit, short) and, for each, -an extra dependency line: -``` -query_data///reads_R1.fastq.gz: genomes/.fna.gz -``` -distinct from the `simulated_data/...` rule for the same specimen. +## Sparse global index -Read count fixed at 100,000 read pairs per genome, independent of genome -size — unlike `simulate_one.sh`'s `simulated_data/` runs, which derive -`n_reads` from a fixed 15x coverage target. A query benchmark does not need -coverage-proportional depth; a fixed pair count keeps the two query runs -comparable to each other and keeps wall/RSS numbers meaningful across -genomes of very different sizes (bacterium vs archaeon). This likely needs a -dedicated `simulate_query_one.sh` (or a parameter to `simulate_one.sh`) -rather than reusing it unchanged, since `n_reads` is currently computed -in-script from genome size. +`global_index_presence_sparse/` is built by `pack_sparse.sh`: copy +`global_index_presence/` wholesale, then `obikmer pack --sparse` in place. +This works directly because `merge`'s pack stage (`merge.rs:252`, +`pack_matrices(false)`) keeps the per-genome column files on disk after +dense-packing — `pack_sparse_bit_matrix` (`obicompactvec/src/bitmatrix/sparse.rs:447`) +reads those, is idempotent, and removes `matrix.pbmx` once the sparse form +is written, so `Persistent::open` falls through to the sparse format +afterward. No separate merge run needed. -**Phase 1 — sparse global index.** New target -`global_index_presence_sparse/index.done`, built from `global_index_presence/` -via `obikmer pack --sparse`. Open question, to verify against the `pack` -implementation before writing the rule: does `pack --sparse` accept an -already dense-packed index in place (`cp -r` + repack), or does it require -the pre-pack column layout, forcing a dedicated merge run instead of reusing -`global_index_presence/`? +## Query runs -**Phase 2 — query runs.** For each of the two `QUERY_SPECIMENS`, run -`obikmer query` against both `global_index_presence` and -`global_index_presence_sparse`, capturing Reporter wall/RSS stats the same -way `merge_presence.sh` does (stderr capture + `parse_reporter`). +`query_one.sh dense|sparse SPECIMEN` runs `obikmer query --count-missing` +against `global_index_presence` or `global_index_presence_sparse`, output +gzipped to `query_{dense,sparse}/SPECIMEN.fasta.gz`, Reporter wall/RSS +captured to `stats/query_{dense,sparse}/SPECIMEN.stats` (same +stderr-parsing convention as `merge_presence.sh`). -**Phase 3 — dense/sparse regression.** `verify_query.py` diffs the two query -JSON outputs per specimen (same matches, same per-genome presence -annotations) → `.stats` CSV (`run,specimen,mismatches,pct`), aggregated by -`aggregate_stats.sh` under a new `query` case. Any mismatch is a real -regression — dense and sparse must be content-identical, only I/O access -pattern differs. +Flags: `--count-missing` only. `--mismatch` is a no-op today +(`query/mod.rs:212-213`, prints "not yet implemented, ignored") — left off +rather than tested for a feature that doesn't exist yet. -**Phase 4 — performance comparison.** No dedicated script: the wall/RSS -columns from Phase 2's `.stats` files, aggregated, are the dense-vs-sparse -performance comparison (the expected win for query on sparse, per the `pack ---sparse` help text). +## Dense/sparse regression -`count` track is out of scope for the sparse branch: `pack --sparse` targets -presence matrices only (per CLI help), no count equivalent confirmed. +`verify_query.py` compares the two query outputs per specimen, matched by +read id (not stream position — the query pipeline chunks input across +worker threads and doesn't guarantee output order). Compares `kmer_count`, +`kmer_missing`, and the full `kmer_strict_matches` map per read. Any +mismatch is a real regression: dense and sparse must be content-identical, +only I/O access pattern differs. `.stats` → `stats/verify_query/`, +aggregated by `aggregate_stats.sh query|verify_query`-style cases +(`query_dense`, `query_sparse`, `verify_query`). + +## Performance comparison + +No dedicated script: the wall/RSS columns from the `query_dense` and +`query_sparse` aggregated `.stats` CSVs are the dense-vs-sparse performance +comparison — the expected win for query on sparse, per the `pack --sparse` +help text. + +## Scope + +`count` track excluded from the sparse branch: `pack --sparse` targets +presence matrices only (per CLI help); `pack_matrices` leaves count +matrices untouched regardless of the `sparse` flag +(`obikindex/src/index.rs:308`). + +## New Makefile targets + +`simulate_query`, `pack_sparse`, `query_dense`, `query_sparse`, +`aggregate_query_dense`, `aggregate_query_sparse`, `verify_query`, +`aggregate_verify_query` — the last three folded into `all`. diff --git a/benchmark/Makefile b/benchmark/Makefile index 5bd9b7a6..be7a1958 100644 --- a/benchmark/Makefile +++ b/benchmark/Makefile @@ -38,6 +38,12 @@ SPECIFIC_PRESENCE_STATS := $(SPECIES:%=stats/specific_kmer_presence/%.stats) SPECIFIC_COUNT_DONE := $(SPECIES:%=specific_index_count/%/index.done) SPECIFIC_COUNT_STATS := $(SPECIES:%=stats/specific_kmer_count/%.stats) SIMULATED_READS := $(foreach s,$(SPECIMENS),simulated_data/$(subst --,/,$s)/reads_R1.fastq.gz) +QUERY_READS := $(foreach s,$(QUERY_SPECIMENS),query_data/$(subst --,/,$s)/reads_R1.fastq.gz) +QUERY_DENSE_DONE := $(QUERY_SPECIMENS:%=query_dense/%.fasta.gz) +QUERY_DENSE_STATS := $(QUERY_SPECIMENS:%=stats/query_dense/%.stats) +QUERY_SPARSE_DONE := $(QUERY_SPECIMENS:%=query_sparse/%.fasta.gz) +QUERY_SPARSE_STATS := $(QUERY_SPECIMENS:%=stats/query_sparse/%.stats) +VERIFY_QUERY_STATS := $(QUERY_SPECIMENS:%=stats/verify_query/%.stats) .NOTPARALLEL: @@ -51,7 +57,11 @@ SIMULATED_READS := $(foreach s,$(SPECIMENS),simulated_data/$(subst --,/,$s)/read aggregate_verify_presence aggregate_verify_count \ verify_merge_presence verify_merge_count \ filter_presence filter_count \ - aggregate_filter_presence aggregate_filter_count + aggregate_filter_presence aggregate_filter_count \ + pack_sparse simulate_query \ + query_dense query_sparse \ + aggregate_query_dense aggregate_query_sparse \ + verify_query aggregate_verify_query verify_merge_presence: stats/verify_merge_presence/current.csv verify_merge_count: stats/verify_merge_count/current.csv @@ -59,12 +69,13 @@ verify_merge_count: stats/verify_merge_count/current.csv all: aggregate_verify_presence aggregate_verify_count \ verify_merge_presence verify_merge_count \ aggregate_filter_presence aggregate_filter_count \ - dist_comparison + dist_comparison \ + aggregate_query_dense aggregate_query_sparse aggregate_verify_query # ── dependency file ─────────────────────────────────────────────────────────── -deps.mk: $(GENOMES) - $(VENV_PY) make_deps.py $^ > $@ +deps.mk: $(GENOMES) make_deps.py + $(VENV_PY) make_deps.py $(GENOMES) > $@ # ── simulation ──────────────────────────────────────────────────────────────── # Prerequisites (genome → reads) are in deps.mk; $< is the genome file. @@ -74,6 +85,14 @@ $(SIMULATED_READS): simulate: $(SIMULATED_READS) +# ── query read simulation (fixed size, independent draw) ─────────────────────── +# Prerequisites (genome → reads) are in deps.mk; $< is the genome file. + +$(QUERY_READS): + bash simulate_query_one.sh $< $(dir $@) + +simulate_query: $(QUERY_READS) + # ── reference kmer sets ─────────────────────────────────────────────────────── # Prerequisites (reads → npz) are in deps.mk. @@ -228,3 +247,40 @@ stats/verify_merge_presence/current.csv: $(REF_NPZS) global_index_presence/index stats/verify_merge_count/current.csv: $(REF_NPZS) global_index_count/index.done bash verify_merge_count.sh + +# ── sparse presence index (query benchmark) ───────────────────────────────────── + +global_index_presence_sparse/index.done: global_index_presence/index.done $(BINARY) + bash pack_sparse.sh + +pack_sparse: global_index_presence_sparse/index.done + +# ── query: dense vs sparse ─────────────────────────────────────────────────────── +# Prerequisites (reads + index → output + .stats) are in deps.mk. + +query_dense/%.fasta.gz \ +stats/query_dense/%.stats &: $(BINARY) + bash query_one.sh dense $* + +query_sparse/%.fasta.gz \ +stats/query_sparse/%.stats &: $(BINARY) + bash query_one.sh sparse $* + +query_dense: $(QUERY_DENSE_DONE) +query_sparse: $(QUERY_SPARSE_DONE) + +aggregate_query_dense: $(QUERY_DENSE_STATS) + bash aggregate_stats.sh query_dense + +aggregate_query_sparse: $(QUERY_SPARSE_STATS) + bash aggregate_stats.sh query_sparse + +# ── query: dense/sparse regression ────────────────────────────────────────────── + +stats/verify_query/%.stats: + bash verify_query_one.sh $* + +verify_query: $(VERIFY_QUERY_STATS) + +aggregate_verify_query: $(VERIFY_QUERY_STATS) + bash aggregate_stats.sh verify_query diff --git a/benchmark/README.md b/benchmark/README.md index 04ad741f..ea683715 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -75,10 +75,41 @@ flowchart TD verify_merge_count --> vmc[("stats/verify_merge_count/")] end + subgraph query ["Query track (2 specimens: E. coli + archaeon)"] + GENOMES --> simulate_query + simulate_query --> qdata[("query_data/")] + + gpres --> pack_sparse + BIN --> pack_sparse + pack_sparse --> gsparse[("global_index_presence_sparse/")] + + qdata --> query_dense + gpres --> query_dense + BIN --> query_dense + query_dense --> qd[("query_dense/")] + query_dense --> qd_stats[("stats/query_dense/")] + qd_stats --> aggregate_query_dense + + qdata --> query_sparse + gsparse --> query_sparse + BIN --> query_sparse + query_sparse --> qs[("query_sparse/")] + query_sparse --> qs_stats[("stats/query_sparse/")] + qs_stats --> aggregate_query_sparse + + qd --> verify_query + qs --> verify_query + verify_query --> vq_stats[("stats/verify_query/")] + vq_stats --> aggregate_verify_query + end + aggregate_verify_presence --> all aggregate_verify_count --> all vmp --> all vmc --> all + aggregate_query_dense --> all + aggregate_query_sparse --> all + aggregate_verify_query --> all all -. "$(MAKE) re-eval" .-> aggregate_filter_presence all -. "$(MAKE) re-eval" .-> aggregate_filter_count ``` @@ -105,6 +136,14 @@ flowchart TD | `aggregate_filter_count` | `aggregate_stats.sh` | Aggregate species-specific kmer stats (count) | | `verify_merge_presence` | `verify_merge_presence.sh` | Verify global presence index against all reference sets | | `verify_merge_count` | `verify_merge_count.sh` | Verify global count index against all reference sets | +| `simulate_query` | `simulate_query_one.sh` | Simulate a fixed-size (100k pairs) read set per query specimen | +| `pack_sparse` | `pack_sparse.sh` | Build `global_index_presence_sparse/` from `global_index_presence/` | +| `query_dense` | `query_one.sh dense` | Query each query specimen's reads against the dense global index | +| `query_sparse` | `query_one.sh sparse` | Query each query specimen's reads against the sparse global index | +| `aggregate_query_dense` | `aggregate_stats.sh` | Aggregate dense query wall/RSS stats | +| `aggregate_query_sparse` | `aggregate_stats.sh` | Aggregate sparse query wall/RSS stats | +| `verify_query` | `verify_query_one.sh` | Diff dense vs sparse query output per specimen (regression check) | +| `aggregate_verify_query` | `aggregate_stats.sh` | Aggregate dense/sparse query regression stats | ## Directory layout @@ -113,13 +152,18 @@ benchmark/ ├── genomes/ # input reference genomes (.fna.gz) ├── simulated_data/ # generated by simulate │ └── // +├── query_data/ # generated by simulate_query (2 specimens, fixed 100k pairs) +│ └── // ├── reference_index/ # reference kmer sets (.npz) ├── specimen_index_presence/ # per-specimen presence indexes ├── specimen_index_count/ # per-specimen count indexes -├── global_index_presence/ # merged global presence index +├── global_index_presence/ # merged global presence index (dense-packed) +├── global_index_presence_sparse/ # global presence index, sparse-packed (query benchmark) ├── global_index_count/ # merged global count index ├── specific_index_presence/ # species-specific presence indexes ├── specific_index_count/ # species-specific count indexes +├── query_dense/ # query output against global_index_presence +├── query_sparse/ # query output against global_index_presence_sparse └── stats/ # all benchmark statistics ├── indexing_presence/ ├── indexing_count/ @@ -128,5 +172,9 @@ benchmark/ ├── specific_kmer_presence/ ├── specific_kmer_count/ ├── verify_merge_presence/ - └── verify_merge_count/ + ├── verify_merge_count/ + ├── pack_sparse/ + ├── query_dense/ + ├── query_sparse/ + └── verify_query/ ``` diff --git a/benchmark/aggregate_stats.sh b/benchmark/aggregate_stats.sh index 19901bb9..89f7fa90 100755 --- a/benchmark/aggregate_stats.sh +++ b/benchmark/aggregate_stats.sh @@ -24,6 +24,12 @@ case "${TYPE}" in specific_kmer_presence|specific_kmer_count) HEADER="run,species,rebuild_wall_s,rebuild_rss_b,pack_wall_s,pack_rss_b,filter_total_wall_s,filter_total_rss_b,select_wall_s,select_rss_b,select_total_wall_s,select_total_rss_b" ;; + query_dense|query_sparse) + HEADER="run,species,strain,query_wall_s,query_rss_b,total_wall_s,total_rss_b" + ;; + verify_query) + HEADER="run,species,strain,n_reads,n_common,missing_in_dense,missing_in_sparse,mismatched,mismatch_pct" + ;; *) echo "ERROR: unknown stats type '${TYPE}'" >&2 exit 1 diff --git a/benchmark/deps.mk b/benchmark/deps.mk index 031dd598..2cf1a49f 100644 --- a/benchmark/deps.mk +++ b/benchmark/deps.mk @@ -197,3 +197,17 @@ specific_index_count/Yersinia_ruckeri/index.done stats/specific_kmer_count/Yersi # Candidozyma_auris specific_index_presence/Candidozyma_auris/index.done stats/specific_kmer_presence/Candidozyma_auris.stats: global_index_presence/index.done specific_index_count/Candidozyma_auris/index.done stats/specific_kmer_count/Candidozyma_auris.stats: global_index_count/index.done + +QUERY_SPECIMENS := Escherichia_coli--K-12_MG1655 Saccharolobus_islandicus--M.16.4 + +# query: Escherichia_coli--K-12_MG1655 +query_data/Escherichia_coli/K-12_MG1655/reads_R1.fastq.gz: genomes/GCF_000005845.2_ASM584v2_genomic.fna.gz +query_dense/Escherichia_coli--K-12_MG1655.fasta.gz stats/query_dense/Escherichia_coli--K-12_MG1655.stats: query_data/Escherichia_coli/K-12_MG1655/reads_R1.fastq.gz global_index_presence/index.done +query_sparse/Escherichia_coli--K-12_MG1655.fasta.gz stats/query_sparse/Escherichia_coli--K-12_MG1655.stats: query_data/Escherichia_coli/K-12_MG1655/reads_R1.fastq.gz global_index_presence_sparse/index.done +stats/verify_query/Escherichia_coli--K-12_MG1655.stats: query_dense/Escherichia_coli--K-12_MG1655.fasta.gz query_sparse/Escherichia_coli--K-12_MG1655.fasta.gz + +# query: Saccharolobus_islandicus--M.16.4 +query_data/Saccharolobus_islandicus/M.16.4/reads_R1.fastq.gz: genomes/GCF_000022445.1_ASM2244v1_genomic.fna.gz +query_dense/Saccharolobus_islandicus--M.16.4.fasta.gz stats/query_dense/Saccharolobus_islandicus--M.16.4.stats: query_data/Saccharolobus_islandicus/M.16.4/reads_R1.fastq.gz global_index_presence/index.done +query_sparse/Saccharolobus_islandicus--M.16.4.fasta.gz stats/query_sparse/Saccharolobus_islandicus--M.16.4.stats: query_data/Saccharolobus_islandicus/M.16.4/reads_R1.fastq.gz global_index_presence_sparse/index.done +stats/verify_query/Saccharolobus_islandicus--M.16.4.stats: query_dense/Saccharolobus_islandicus--M.16.4.fasta.gz query_sparse/Saccharolobus_islandicus--M.16.4.fasta.gz diff --git a/benchmark/global_index_presence_sparse/.obikmer.lock b/benchmark/global_index_presence_sparse/.obikmer.lock new file mode 100644 index 00000000..e69de29b diff --git a/benchmark/make_deps.py b/benchmark/make_deps.py index 03f7e2a3..a0381a70 100644 --- a/benchmark/make_deps.py +++ b/benchmark/make_deps.py @@ -13,6 +13,11 @@ STOP_WORDS = {'complete', 'chromosome', 'whole', 'sequence', 'genome', 'endosymbiont', 'of'} STOP_PREFIXES = ('scaffold', 'contig', 'plasmid') +# Specimens used as read sources for the query benchmark (see +# DevDocMD/implementation/benchmark_query_testing.md): one common bacterium, +# one distant lineage (the only archaeon in SPECIES). +QUERY_SPECIMENS = ['Escherichia_coli--K-12_MG1655', 'Saccharolobus_islandicus--M.16.4'] + def is_stop(tok): t = tok.lower() @@ -113,6 +118,27 @@ def main(): print(f'{sp_done} {sp_stats}: global_index_presence/index.done') print(f'{sc_done} {sc_stats}: global_index_count/index.done') + print() + print('QUERY_SPECIMENS :=', ' '.join(QUERY_SPECIMENS)) + + by_specimen = {e[0]: e for e in entries} + for specimen in QUERY_SPECIMENS: + _, species, sim_dir, genome = by_specimen[specimen] + query_dir = sim_dir.replace('simulated_data/', 'query_data/', 1) + reads = f'{query_dir}/reads_R1.fastq.gz' + dense_out = f'query_dense/{specimen}.fasta.gz' + dense_stat = f'stats/query_dense/{specimen}.stats' + sparse_out = f'query_sparse/{specimen}.fasta.gz' + sparse_stat = f'stats/query_sparse/{specimen}.stats' + vq_stat = f'stats/verify_query/{specimen}.stats' + + print() + print(f'# query: {specimen}') + print(f'{reads}: {genome}') + print(f'{dense_out} {dense_stat}: {reads} global_index_presence/index.done') + print(f'{sparse_out} {sparse_stat}: {reads} global_index_presence_sparse/index.done') + print(f'{vq_stat}: {dense_out} {sparse_out}') + if __name__ == '__main__': main() diff --git a/benchmark/pack_sparse.sh b/benchmark/pack_sparse.sh new file mode 100755 index 00000000..47ea2416 --- /dev/null +++ b/benchmark/pack_sparse.sh @@ -0,0 +1,86 @@ +#!/usr/bin/env bash +# Builds global_index_presence_sparse/ from global_index_presence/ by +# copying the index (column files are kept in place after merge's dense +# pack — see obikindex::KmerIndex::pack_matrices) and repacking in place +# with --sparse. +# Outputs: +# global_index_presence_sparse/index.done (copied from source) +# stats/pack_sparse/current.stats (one CSV data row, no header) +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +BINARY="${SCRIPT_DIR}/../src/target/release/obikmer" +SOURCE="${SCRIPT_DIR}/global_index_presence" +OUTPUT="${SCRIPT_DIR}/global_index_presence_sparse" +STATS_DIR="${SCRIPT_DIR}/stats/pack_sparse" +STATS_FILE="${STATS_DIR}/current.stats" + +mkdir -p "${STATS_DIR}" + +echo "[pack_sparse] ${SOURCE} → ${OUTPUT}" + +rm -rf "${OUTPUT}" +cp -r "${SOURCE}" "${OUTPUT}" + +STDERR_LOG=$(mktemp) +trap 'rm -f "${STDERR_LOG}"' EXIT + +"${BINARY}" pack --sparse "${OUTPUT}" 2>"${STDERR_LOG}" + +cat "${STDERR_LOG}" >&2 + +python3 - "${STDERR_LOG}" <<'PYEOF' >"${STATS_FILE}" +import sys, re + +logfile = sys.argv[1] + +def strip_ansi(s): + return re.sub(r'\x1b\[[\x30-\x3f]*[\x20-\x2f]*[\x40-\x7e]', '', s) + +def parse_wall(s): + s = s.strip() + if s.endswith('ms'): return float(s[:-2]) / 1000.0 + if s.endswith('s'): return float(s[:-1]) + return 0.0 + +def parse_rss(s): + m = re.match(r'([\d.]+)\s*(GB|MB|KB|B)', s.strip()) + if not m: return 0 + return int(float(m.group(1)) * {'GB': 1<<30, 'MB': 1<<20, 'KB': 1024, 'B': 1}[m.group(2)]) + +def is_sep(s): + return bool(s) and not re.search(r'[A-Za-z0-9]', s) + +stats = {} +state = 'scan' +with open(logfile, errors='replace') as fh: + for raw in fh: + line = strip_ansi(raw.rstrip('\n')) + s = line.strip() + if state == 'scan': + if re.search(r'\bstage\b.*\bwall\b', line): + state = 'in_header' + elif state == 'in_header': + if is_sep(s): state = 'rows' + elif state == 'rows': + if is_sep(s): state = 'total' + elif s: + parts = re.split(r' +', s) + if len(parts) >= 4: + stats[parts[0]] = (parse_wall(parts[1]), parse_rss(parts[3])) + elif state == 'total': + if s: + parts = re.split(r' +', s) + if len(parts) >= 3: + stats['TOTAL'] = (parse_wall(parts[1]), + parse_rss(parts[3]) if len(parts) > 3 else 0) + break + +w, r = stats.get('pack', ('', '')) +tw, tr = stats.get('TOTAL', ('', '')) +row = [f'{w:.3f}' if isinstance(w, float) else '', str(r), + f'{tw:.3f}' if isinstance(tw, float) else '', str(tr)] +print(','.join(row)) +PYEOF + +echo "Done → ${OUTPUT}" diff --git a/benchmark/query_one.sh b/benchmark/query_one.sh new file mode 100755 index 00000000..a984c7e1 --- /dev/null +++ b/benchmark/query_one.sh @@ -0,0 +1,105 @@ +#!/usr/bin/env bash +# Usage: query_one.sh MODE SPECIMEN +# MODE = dense | sparse +# SPECIMEN = "species--strain" (Make pattern stem), reads from query_data/ +# Outputs: +# query_MODE/SPECIMEN.fasta.gz (obikmer query output, --count-missing) +# stats/query_MODE/SPECIMEN.stats (one CSV data row, no header) +set -euo pipefail + +MODE="$1" +SPECIMEN="$2" +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +BINARY="${SCRIPT_DIR}/../src/target/release/obikmer" + +case "${MODE}" in + dense) INDEX="${SCRIPT_DIR}/global_index_presence" ;; + sparse) INDEX="${SCRIPT_DIR}/global_index_presence_sparse" ;; + *) echo "ERROR: unknown mode '${MODE}' (expected dense|sparse)" >&2; exit 1 ;; +esac + +species="${SPECIMEN%%--*}" +strain="${SPECIMEN#*--}" + +READS_DIR="${SCRIPT_DIR}/query_data/${species}/${strain}" +OUT_DIR="${SCRIPT_DIR}/query_${MODE}" +STATS_DIR="${SCRIPT_DIR}/stats/query_${MODE}" +OUT_FILE="${OUT_DIR}/${SPECIMEN}.fasta.gz" +STATS_FILE="${STATS_DIR}/${SPECIMEN}.stats" + +mkdir -p "${OUT_DIR}" "${STATS_DIR}" + +r1="${READS_DIR}/reads_R1.fastq.gz" +r2="${READS_DIR}/reads_R2.fastq.gz" +if [[ ! -f "${r1}" || ! -f "${r2}" ]]; then + echo "ERROR: reads not found in ${READS_DIR}" >&2 + exit 1 +fi + +echo "[${SPECIMEN}] query (${MODE}) → ${OUT_FILE}" + +STDERR_LOG=$(mktemp) +trap 'rm -f "${STDERR_LOG}"' EXIT + +"${BINARY}" query \ + --count-missing \ + "${INDEX}" "${r1}" "${r2}" \ + 2>"${STDERR_LOG}" \ + | gzip >"${OUT_FILE}" + +cat "${STDERR_LOG}" >&2 + +python3 - "${species}" "${strain}" "${STDERR_LOG}" <<'PYEOF' >"${STATS_FILE}" +import sys, re + +species, strain, logfile = sys.argv[1], sys.argv[2], sys.argv[3] + +def strip_ansi(s): + return re.sub(r'\x1b\[[\x30-\x3f]*[\x20-\x2f]*[\x40-\x7e]', '', s) + +def parse_wall(s): + s = s.strip() + if s.endswith('ms'): return float(s[:-2]) / 1000.0 + if s.endswith('s'): return float(s[:-1]) + return 0.0 + +def parse_rss(s): + m = re.match(r'([\d.]+)\s*(GB|MB|KB|B)', s.strip()) + if not m: return 0 + return int(float(m.group(1)) * {'GB': 1<<30, 'MB': 1<<20, 'KB': 1024, 'B': 1}[m.group(2)]) + +def is_sep(s): + return bool(s) and not re.search(r'[A-Za-z0-9]', s) + +stats = {} +state = 'scan' +with open(logfile, errors='replace') as fh: + for raw in fh: + line = strip_ansi(raw.rstrip('\n')) + s = line.strip() + if state == 'scan': + if re.search(r'\bstage\b.*\bwall\b', line): + state = 'in_header' + elif state == 'in_header': + if is_sep(s): state = 'rows' + elif state == 'rows': + if is_sep(s): state = 'total' + elif s: + parts = re.split(r' +', s) + if len(parts) >= 4: + stats[parts[0]] = (parse_wall(parts[1]), parse_rss(parts[3])) + elif state == 'total': + if s: + parts = re.split(r' +', s) + if len(parts) >= 3: + stats['TOTAL'] = (parse_wall(parts[1]), + parse_rss(parts[3]) if len(parts) > 3 else 0) + break + +qw, qr = stats.get('query', ('', '')) +tw, tr = stats.get('TOTAL', ('', '')) +row = [species, strain, + f'{qw:.3f}' if isinstance(qw, float) else '', str(qr), + f'{tw:.3f}' if isinstance(tw, float) else '', str(tr)] +print(','.join(row)) +PYEOF diff --git a/benchmark/simulate_query_one.sh b/benchmark/simulate_query_one.sh new file mode 100755 index 00000000..dd567347 --- /dev/null +++ b/benchmark/simulate_query_one.sh @@ -0,0 +1,34 @@ +#!/usr/bin/env bash +# Usage: simulate_query_one.sh genome.fna.gz output_dir +# Simulates a fixed-size paired-end HiSeq read set for the query benchmark. +# Unlike simulate_one.sh (coverage-proportional, used to build the indexed +# specimens), this always draws N_READS pairs regardless of genome size — +# query benchmark numbers (wall/RSS) must stay comparable across genomes of +# very different sizes. Independent iss run (unseeded), so error draw +# differs from any simulated_data/ reads for the same genome. +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +ISS="${SCRIPT_DIR}/../.venv/bin/iss" +N_READS=100000 +CPUS="${CPUS:-$(sysctl -n hw.logicalcpu 2>/dev/null || nproc 2>/dev/null || echo 2)}" + +genome_file="$1" +out_dir="$2" + +mkdir -p "${out_dir}" + +tmp_fasta=$(mktemp "${TMPDIR:-/tmp}/obikmer_XXXXXX.fna") +trap 'rm -f "${tmp_fasta}"' EXIT + +gzip -dc "${genome_file}" > "${tmp_fasta}" + +echo "[${out_dir}] ${N_READS} read pairs (query benchmark, fixed size)" + +"${ISS}" generate \ + --genomes "${tmp_fasta}" \ + --model HiSeq \ + --n_reads "${N_READS}" \ + --cpus "${CPUS}" \ + --compress \ + --output "${out_dir}/reads" diff --git a/benchmark/verify_query.py b/benchmark/verify_query.py new file mode 100755 index 00000000..cc7e8a11 --- /dev/null +++ b/benchmark/verify_query.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python3 +"""Compare dense vs sparse obikmer query output for one specimen. + +Both files are `obikmer query --count-missing` output: OBITools4-style +FASTA, one record per read — `>id {"kmer_count":N,"kmer_missing":M, +"kmer_strict_matches":{"label":count,...}}`. Packing format (dense vs +sparse presence matrix) must not change query results — only I/O access +pattern differs. Matched by read id rather than by stream position: the +query pipeline processes input in chunks across worker threads and does +not guarantee output order matches input order. + +Output to stdout: one CSV row + species, strain, n_reads, n_common, missing_in_dense, missing_in_sparse, mismatched, mismatch_pct +""" +import argparse +import gzip +import json +import sys + + +def load(path): + records = {} + opener = gzip.open if path.endswith('.gz') else open + with opener(path, 'rt') as fh: + for line in fh: + if not line.startswith('>'): + continue + header = line[1:].rstrip('\n') + read_id, _, json_part = header.partition(' ') + records[read_id] = json.loads(json_part) if json_part else {} + return records + + +def annotations_equal(a, b): + return ( + a.get('kmer_count') == b.get('kmer_count') + and a.get('kmer_missing') == b.get('kmer_missing') + and a.get('kmer_strict_matches', {}) == b.get('kmer_strict_matches', {}) + ) + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument('--species', required=True) + ap.add_argument('--strain', required=True) + ap.add_argument('dense_fasta') + ap.add_argument('sparse_fasta') + args = ap.parse_args() + + dense = load(args.dense_fasta) + sparse = load(args.sparse_fasta) + + dense_ids, sparse_ids = set(dense), set(sparse) + common = dense_ids & sparse_ids + missing_in_dense = len(sparse_ids - dense_ids) + missing_in_sparse = len(dense_ids - sparse_ids) + + mismatched = sum(1 for rid in common if not annotations_equal(dense[rid], sparse[rid])) + + n_reads = len(dense_ids | sparse_ids) + mismatch_pct = 100.0 * (mismatched + missing_in_dense + missing_in_sparse) / n_reads if n_reads else 0.0 + + print(f'{args.species},{args.strain},{n_reads},{len(common)},' + f'{missing_in_dense},{missing_in_sparse},{mismatched},{mismatch_pct:.6f}') + + +if __name__ == '__main__': + main() diff --git a/benchmark/verify_query_one.sh b/benchmark/verify_query_one.sh new file mode 100755 index 00000000..b416e634 --- /dev/null +++ b/benchmark/verify_query_one.sh @@ -0,0 +1,28 @@ +#!/usr/bin/env bash +# Usage: verify_query_one.sh SPECIMEN +# SPECIMEN = "species--strain" (Make pattern stem) +# Output: stats/verify_query/SPECIMEN.stats (one CSV data row, no header) +set -euo pipefail + +SPECIMEN="$1" +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +PYTHON="${SCRIPT_DIR}/../.venv/bin/python3" +VERIFY_PY="${SCRIPT_DIR}/verify_query.py" + +species="${SPECIMEN%%--*}" +strain="${SPECIMEN#*--}" + +DENSE="${SCRIPT_DIR}/query_dense/${SPECIMEN}.fasta.gz" +SPARSE="${SCRIPT_DIR}/query_sparse/${SPECIMEN}.fasta.gz" +STATS_DIR="${SCRIPT_DIR}/stats/verify_query" +STATS_FILE="${STATS_DIR}/${SPECIMEN}.stats" + +mkdir -p "${STATS_DIR}" + +echo "[${SPECIMEN}] verifying query (dense vs sparse)" + +"${PYTHON}" "${VERIFY_PY}" \ + --species "${species}" \ + --strain "${strain}" \ + "${DENSE}" "${SPARSE}" \ + >"${STATS_FILE}" diff --git a/mkdocs.yml b/mkdocs.yml index 0731ce69..b1660442 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -55,7 +55,7 @@ nav: - Kmer filtering: implementation/filtering.md - Select command: implementation/select.md - obitaxonomy crate: implementation/obitaxonomy.md - - "Benchmark: query-path testing (discussion)": implementation/benchmark_query_testing.md + - "Benchmark: query-path testing": implementation/benchmark_query_testing.md - Architecture: - Sequences: architecture/sequences/invariant.md - Kmer index: architecture/index_architecture.md