feat: add batched int group stats API and expand benchmark variants

Introduces a `batch_int_group_stats` API for computing presence counts, sums, minimums, and maximums across sparse and dense matrix representations. The selection layer now utilizes this batched approach to optimize aggregation semantics for boolean and numeric operations. Additionally, reorganizes the benchmarking infrastructure to support querying across presence and count index variants in both dense and sparse formats, including new packing scripts and updated statistics aggregation.
This commit is contained in:
Eric Coissac
2026-08-28 23:37:16 +02:00
parent 579cfd8752
commit 54e4600120
13 changed files with 546 additions and 163 deletions
+68 -34
View File
@@ -39,10 +39,14 @@ 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)
QUERY_PRESENCE_DENSE_DONE := $(QUERY_SPECIMENS:%=query_presence_dense/%.fasta.gz)
QUERY_PRESENCE_DENSE_STATS := $(QUERY_SPECIMENS:%=stats/query_presence_dense/%.stats)
QUERY_PRESENCE_SPARSE_DONE := $(QUERY_SPECIMENS:%=query_presence_sparse/%.fasta.gz)
QUERY_PRESENCE_SPARSE_STATS := $(QUERY_SPECIMENS:%=stats/query_presence_sparse/%.stats)
QUERY_COUNT_DENSE_DONE := $(QUERY_SPECIMENS:%=query_count_dense/%.fasta.gz)
QUERY_COUNT_DENSE_STATS := $(QUERY_SPECIMENS:%=stats/query_count_dense/%.stats)
QUERY_COUNT_SPARSE_DONE := $(QUERY_SPECIMENS:%=query_count_sparse/%.fasta.gz)
QUERY_COUNT_SPARSE_STATS := $(QUERY_SPECIMENS:%=stats/query_count_sparse/%.stats)
VERIFY_QUERY_STATS := $(QUERY_SPECIMENS:%=stats/verify_query/%.stats)
.NOTPARALLEL:
@@ -58,9 +62,11 @@ VERIFY_QUERY_STATS := $(QUERY_SPECIMENS:%=stats/verify_query/%.stats)
verify_merge_presence verify_merge_count \
filter_presence filter_count \
aggregate_filter_presence aggregate_filter_count \
pack_sparse simulate_query \
query_dense query_sparse \
aggregate_query_dense aggregate_query_sparse \
pack_dense_presence pack_dense_count simulate_query \
query_presence_dense query_presence_sparse \
query_count_dense query_count_sparse \
aggregate_query_presence_dense aggregate_query_presence_sparse \
aggregate_query_count_dense aggregate_query_count_sparse \
verify_query aggregate_verify_query
verify_merge_presence: stats/verify_merge_presence/current.csv
@@ -70,7 +76,9 @@ all: aggregate_verify_presence aggregate_verify_count \
verify_merge_presence verify_merge_count \
aggregate_filter_presence aggregate_filter_count \
dist_comparison \
aggregate_query_dense aggregate_query_sparse aggregate_verify_query
aggregate_query_presence_dense aggregate_query_presence_sparse \
aggregate_query_count_dense aggregate_query_count_sparse \
aggregate_verify_query
# ── dependency file ───────────────────────────────────────────────────────────
@@ -114,11 +122,11 @@ $(OBIKMER_PRESENCE_DIST) &: global_index_presence/index.done $(BINARY)
mkdir -p obikmer_dist/presence
$(BINARY) phylo \
--output obikmer_dist/presence/jaccard \
--metric jaccard --shared-kmers --nj \
--distance jaccard --csv --shared-kmers --nj \
global_index_presence
$(BINARY) phylo \
--output obikmer_dist/presence/hamming \
--metric hamming --nj \
--distance hamming --csv --nj \
global_index_presence
obikmer_dist_presence: $(OBIKMER_PRESENCE_DIST)
@@ -129,31 +137,31 @@ $(OBIKMER_COUNT_DIST) &: global_index_count/index.done $(BINARY)
mkdir -p obikmer_dist/count
$(BINARY) phylo \
--output obikmer_dist/count/jaccard \
--metric jaccard --shared-kmers --nj \
--distance jaccard --csv --shared-kmers --nj \
global_index_count
$(BINARY) phylo \
--output obikmer_dist/count/bray_curtis \
--metric bray-curtis --nj \
--distance bray-curtis --csv --nj \
global_index_count
$(BINARY) phylo \
--output obikmer_dist/count/relfreq_bray_curtis \
--metric relfreq-bray-curtis --nj \
--distance relfreq-bray-curtis --csv --nj \
global_index_count
$(BINARY) phylo \
--output obikmer_dist/count/euclidean \
--metric euclidean --nj \
--distance euclidean --csv --nj \
global_index_count
$(BINARY) phylo \
--output obikmer_dist/count/relfreq_euclidean \
--metric relfreq-euclidean --nj \
--distance relfreq-euclidean --csv --nj \
global_index_count
$(BINARY) phylo \
--output obikmer_dist/count/hellinger \
--metric hellinger --nj \
--distance hellinger --csv --nj \
global_index_count
$(BINARY) phylo \
--output obikmer_dist/count/hellinger_euclidean \
--metric hellinger-euclidean --nj \
--distance hellinger-euclidean --csv --nj \
global_index_count
obikmer_dist_count: $(OBIKMER_COUNT_DIST)
@@ -248,32 +256,58 @@ 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) ─────────────────────────────────────
# ── dense variants (query benchmark) ────────────────────────────────────────────
# `merge` packs sparse by default (2026-08-28) — global_index_presence/
# global_index_count *are* the sparse variants already; the dense ones are
# built explicitly here, from a hard-link-based copy (see
# copy_index_hardlink.sh) rather than a full `cp -r`.
global_index_presence_sparse/index.done: global_index_presence/index.done $(BINARY)
bash pack_sparse.sh
global_index_presence_dense/index.done: global_index_presence/index.done $(BINARY)
bash pack_dense.sh presence
pack_sparse: global_index_presence_sparse/index.done
# Rebuilt from the per-specimen count sources directly (via `merge --dense`),
# not repacked from global_index_count — see pack_dense.sh's own comment.
global_index_count_dense/index.done: $(COUNT_DONE) $(BINARY)
bash pack_dense.sh count
# ── query: dense vs sparse ───────────────────────────────────────────────────────
pack_dense_presence: global_index_presence_dense/index.done
pack_dense_count: global_index_count_dense/index.done
# ── query: dense vs sparse, presence and count ──────────────────────────────────
# Prerequisites (reads + index → output + .stats) are in deps.mk.
query_dense/%.fasta.gz \
stats/query_dense/%.stats &: $(BINARY)
bash query_one.sh dense $*
query_presence_dense/%.fasta.gz \
stats/query_presence_dense/%.stats &: $(BINARY) global_index_presence_dense/index.done
bash query_one.sh presence dense $*
query_sparse/%.fasta.gz \
stats/query_sparse/%.stats &: $(BINARY)
bash query_one.sh sparse $*
query_presence_sparse/%.fasta.gz \
stats/query_presence_sparse/%.stats &: $(BINARY) global_index_presence/index.done
bash query_one.sh presence sparse $*
query_dense: $(QUERY_DENSE_DONE)
query_sparse: $(QUERY_SPARSE_DONE)
query_count_dense/%.fasta.gz \
stats/query_count_dense/%.stats &: $(BINARY) global_index_count_dense/index.done
bash query_one.sh count dense $*
aggregate_query_dense: $(QUERY_DENSE_STATS)
bash aggregate_stats.sh query_dense
query_count_sparse/%.fasta.gz \
stats/query_count_sparse/%.stats &: $(BINARY) global_index_count/index.done
bash query_one.sh count sparse $*
aggregate_query_sparse: $(QUERY_SPARSE_STATS)
bash aggregate_stats.sh query_sparse
query_presence_dense: $(QUERY_PRESENCE_DENSE_DONE)
query_presence_sparse: $(QUERY_PRESENCE_SPARSE_DONE)
query_count_dense: $(QUERY_COUNT_DENSE_DONE)
query_count_sparse: $(QUERY_COUNT_SPARSE_DONE)
aggregate_query_presence_dense: $(QUERY_PRESENCE_DENSE_STATS)
bash aggregate_stats.sh query_presence_dense
aggregate_query_presence_sparse: $(QUERY_PRESENCE_SPARSE_STATS)
bash aggregate_stats.sh query_presence_sparse
aggregate_query_count_dense: $(QUERY_COUNT_DENSE_STATS)
bash aggregate_stats.sh query_count_dense
aggregate_query_count_sparse: $(QUERY_COUNT_SPARSE_STATS)
bash aggregate_stats.sh query_count_sparse
# ── query: dense/sparse regression ──────────────────────────────────────────────
+4 -2
View File
@@ -1,6 +1,8 @@
#!/usr/bin/env bash
# Usage: aggregate_stats.sh TYPE
# TYPE = indexing_presence | indexing_count | verify_presence | verify_count
# TYPE = indexing_presence | indexing_count | verify_presence | verify_count |
# query_presence_dense | query_presence_sparse |
# query_count_dense | query_count_sparse
#
# Reads all stats/TYPE/*.stats files (one CSV data row each, no header).
# Creates a new stats/TYPE/run_NNN.csv only if any .stats file is newer than
@@ -24,7 +26,7 @@ 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)
query_presence_dense|query_presence_sparse|query_count_dense|query_count_sparse)
HEADER="run,species,strain,query_wall_s,query_rss_b,total_wall_s,total_rss_b"
;;
verify_query)
+116
View File
@@ -0,0 +1,116 @@
#!/usr/bin/env bash
# Usage: pack_dense.sh KIND (KIND = presence | count)
#
# Builds global_index_KIND_dense/ from global_index_KIND/ — `merge` packs
# sparse by default (2026-08-28), so the dense comparison arm needs an
# explicit rebuild.
#
# `obikmer pack --dense` cannot do this: it only converts the raw, freshly-
# built columnar (per-genome-file, unpacked) matrix into a packed format —
# `finalize_indexed` already packs (sparse by default) as the last step of
# every index-building command, so there is no columnar leftover for a
# second `pack` invocation to work from; it fails ("No such file or
# directory", `obicompactvec::bitmatrix::packed::pack_bit_matrix` looking
# for a `meta.json` that packing already cleaned up).
#
# `obikmer select` doesn't have that limitation — it always rebuilds its
# output from scratch via `MatrixBuilder`, reading the source through the
# format-agnostic `PersistentBitMatrix`/`PersistentIntMatrix` (Sparse
# included, both content kinds — see `obicompactvec::batch_presence_counts`/
# `batch_int_group_stats`) — so a full, unaggregated passthrough (`--select`
# naming every genome, no `--group`) with `--dense` genuinely repacks
# Sparse → Dense for either kind, and gets the layer-identity hard-linking
# already implemented in `obikselect::select_layer::copy_layer_files` for
# free.
#
# Outputs:
# global_index_KIND_dense/index.done (rebuilt via `select`, dense-packed)
# stats/pack_dense_KIND/current.stats (one CSV data row, no header)
set -euo pipefail
KIND="$1"
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
BINARY="${SCRIPT_DIR}/../src/target/release/obikmer"
SOURCE="${SCRIPT_DIR}/global_index_${KIND}"
OUTPUT="${SCRIPT_DIR}/global_index_${KIND}_dense"
STATS_DIR="${SCRIPT_DIR}/stats/pack_dense_${KIND}"
STATS_FILE="${STATS_DIR}/current.stats"
mkdir -p "${STATS_DIR}"
echo "[pack_dense_${KIND}] ${SOURCE}${OUTPUT}"
LABELS=$("${BINARY}" annotate "${SOURCE}" --dump | tail -n +2 | python3 -c "
import sys, csv
r = csv.reader(sys.stdin)
print(','.join(row[0] for row in r if row))
")
STDERR_LOG=$(mktemp)
trap 'rm -f "${STDERR_LOG}"' EXIT
"${BINARY}" select \
--output "${OUTPUT}" \
--force \
--dense \
--select "${LABELS}" \
"${SOURCE}" \
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('select', ('', ''))
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}"
-86
View File
@@ -1,86 +0,0 @@
#!/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}"
+20 -11
View File
@@ -1,20 +1,29 @@
#!/usr/bin/env bash
# Usage: query_one.sh MODE SPECIMEN
# MODE = dense | sparse
# Usage: query_one.sh KIND MODE SPECIMEN
# KIND = presence | count
# MODE = dense | sparse — `merge` packs sparse by default (2026-08-28), so
# "sparse" is global_index_KIND itself, unmodified; "dense" is the
# explicitly repacked global_index_KIND_dense (see pack_dense.sh).
# 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)
# query_KIND_MODE/SPECIMEN.fasta.gz (obikmer query output, --count-missing)
# stats/query_KIND_MODE/SPECIMEN.stats (one CSV data row, no header)
set -euo pipefail
MODE="$1"
SPECIMEN="$2"
KIND="$1"
MODE="$2"
SPECIMEN="$3"
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
BINARY="${SCRIPT_DIR}/../src/target/release/obikmer"
case "${KIND}" in
presence|count) ;;
*) echo "ERROR: unknown kind '${KIND}' (expected presence|count)" >&2; exit 1 ;;
esac
case "${MODE}" in
dense) INDEX="${SCRIPT_DIR}/global_index_presence" ;;
sparse) INDEX="${SCRIPT_DIR}/global_index_presence_sparse" ;;
sparse) INDEX="${SCRIPT_DIR}/global_index_${KIND}" ;;
dense) INDEX="${SCRIPT_DIR}/global_index_${KIND}_dense" ;;
*) echo "ERROR: unknown mode '${MODE}' (expected dense|sparse)" >&2; exit 1 ;;
esac
@@ -22,8 +31,8 @@ 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_DIR="${SCRIPT_DIR}/query_${KIND}_${MODE}"
STATS_DIR="${SCRIPT_DIR}/stats/query_${KIND}_${MODE}"
OUT_FILE="${OUT_DIR}/${SPECIMEN}.fasta.gz"
STATS_FILE="${STATS_DIR}/${SPECIMEN}.stats"
@@ -36,7 +45,7 @@ if [[ ! -f "${r1}" || ! -f "${r2}" ]]; then
exit 1
fi
echo "[${SPECIMEN}] query (${MODE}) → ${OUT_FILE}"
echo "[${SPECIMEN}] query (${KIND}/${MODE}) → ${OUT_FILE}"
STDERR_LOG=$(mktemp)
trap 'rm -f "${STDERR_LOG}"' EXIT
+2 -2
View File
@@ -12,8 +12,8 @@ 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"
DENSE="${SCRIPT_DIR}/query_presence_dense/${SPECIMEN}.fasta.gz"
SPARSE="${SCRIPT_DIR}/query_presence_sparse/${SPECIMEN}.fasta.gz"
STATS_DIR="${SCRIPT_DIR}/stats/verify_query"
STATS_FILE="${STATS_DIR}/${SPECIMEN}.stats"