diff --git a/.gitignore b/.gitignore
index 7337b94c..c99af57f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -34,12 +34,14 @@ benchmark/genomes
 benchmark/genomes_orig
 benchmark/simulated_data
 benchmark/specimen_index_presence
-benchmark/specimen_index_count
+benchmark/global_index_count_dense
 benchmark/global_index_presence
 benchmark/global_index_presence_orig
 benchmark/global_index_presence_sav
+benchmark/global_index_presence_dense
 benchmark/all_specific
 benchmark/global_index_count
+benchmark/specimen_index_count
 benchmark/stats
 benchmark/reference_index
 benchmark/reference_dist
diff --git a/benchmark/Makefile b/benchmark/Makefile
index be7a1958..b149e41d 100644
--- a/benchmark/Makefile
+++ b/benchmark/Makefile
@@ -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 ──────────────────────────────────────────────
 
diff --git a/benchmark/aggregate_stats.sh b/benchmark/aggregate_stats.sh
index 89f7fa90..2e668be1 100755
--- a/benchmark/aggregate_stats.sh
+++ b/benchmark/aggregate_stats.sh
@@ -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)
diff --git a/benchmark/pack_dense.sh b/benchmark/pack_dense.sh
new file mode 100755
index 00000000..0d747c4a
--- /dev/null
+++ b/benchmark/pack_dense.sh
@@ -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}"
diff --git a/benchmark/pack_sparse.sh b/benchmark/pack_sparse.sh
deleted file mode 100755
index 47ea2416..00000000
--- a/benchmark/pack_sparse.sh
+++ /dev/null
@@ -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}"
diff --git a/benchmark/query_one.sh b/benchmark/query_one.sh
index a984c7e1..093ff8e6 100755
--- a/benchmark/query_one.sh
+++ b/benchmark/query_one.sh
@@ -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
diff --git a/benchmark/verify_query_one.sh b/benchmark/verify_query_one.sh
index b416e634..ba00dc24 100755
--- a/benchmark/verify_query_one.sh
+++ b/benchmark/verify_query_one.sh
@@ -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"
 
diff --git a/src/obicompactvec/src/int_group_ops.rs b/src/obicompactvec/src/int_group_ops.rs
new file mode 100644
index 00000000..d0091da1
--- /dev/null
+++ b/src/obicompactvec/src/int_group_ops.rs
@@ -0,0 +1,163 @@
+//! Batched, row-major-friendly group statistics for [`PersistentIntMatrix`]
+//! — the count-matrix counterpart of
+//! [`crate::bitmatrix::batch_presence_counts`]. Exists for the same two
+//! reasons: `Sparse` has no on-disk column representation (`col_view` is a
+//! deliberate `panic!` there), so a real column can only be synthesised by
+//! scanning every row — exactly what this does directly, via
+//! [`PersistentSparseCompactIntMatrix::for_each_cell_in_row`] — and even
+//! for `Columnar`/`Packed`, computing every group's stats together reads a
+//! column shared by several groups once, not once per group.
+//!
+//! Unlike the bit-matrix case, a count matrix's `sum`/`min`/`max` are
+//! genuine per-value reductions, not derivable from a single presence
+//! count — so this tracks four running quantities per group instead of
+//! one: `presence_count` (cells `>= threshold`), `sum`, `min`, `max`.
+//! `min`/`max` need one more piece of care a dense `col_view` never has to:
+//! a column absent from a sparse row is implicitly `0`, exactly like a
+//! `col_view` row that was never written to — so `max` (0 is never a new
+//! maximum once any real value has been seen) needs no correction, but
+//! `min` does: if fewer than the full group was actually present at a row,
+//! at least one member was implicitly `0`, so the true minimum is `0`
+//! regardless of what was seen among the present ones. Tracked via a fifth,
+//! internal-only accumulator (`present_count`, distinct from
+//! `presence_count`: the former is "how many group members exist at this
+//! row at all", the latter is "how many clear the threshold") and applied
+//! as a single O(n) correction pass per group once the main scan is done.
+//!
+//! `threshold == 0` is `--presence-threshold`'s CLI default, and needs its
+//! own case: every cell, present or implicitly absent, satisfies `v >= 0`,
+//! so `presence_count` is trivially `group.len()` at every row — computed
+//! directly, without paying for the scan's threshold check at all (see
+//! `PersistentIntMatrix::partial_group_presence_count`'s own `col_view`
+//! based implementation, which reaches the same result by iterating a
+//! dense view that materialises the implicit zeros explicitly).
+
+use std::collections::HashMap;
+use std::io;
+
+use crate::colgroup::ColGroup;
+use crate::intmatrix::PersistentIntMatrix;
+use crate::tempintvec::{TempCompactIntVec, TempCompactIntVecBuilder};
+
+/// One group's batched statistics — see the module docs for exactly what
+/// each field means and how `min` differs from a plain per-cell minimum.
+pub struct IntGroupStats {
+    pub presence_count: TempCompactIntVec,
+    pub sum: TempCompactIntVec,
+    pub min: TempCompactIntVec,
+    pub max: TempCompactIntVec,
+}
+
+struct GroupAccum {
+    present_count: TempCompactIntVecBuilder,
+    presence_count: TempCompactIntVecBuilder,
+    sum: TempCompactIntVecBuilder,
+    min: TempCompactIntVecBuilder,
+    max: TempCompactIntVecBuilder,
+}
+
+impl GroupAccum {
+    fn new(n: usize) -> io::Result<Self> {
+        Ok(Self {
+            present_count: TempCompactIntVecBuilder::new(n)?,
+            presence_count: TempCompactIntVecBuilder::new(n)?,
+            sum: TempCompactIntVecBuilder::new(n)?,
+            min: TempCompactIntVecBuilder::new(n)?,
+            max: TempCompactIntVecBuilder::new(n)?,
+        })
+    }
+
+    /// One cell (`value` at `slot`) belonging to this group.
+    #[inline]
+    fn touch(&mut self, slot: usize, value: u32, threshold: u32) {
+        let seen_before = self.present_count.get(slot);
+        self.present_count.set(slot, seen_before + 1);
+        if value >= threshold {
+            self.presence_count.set(slot, self.presence_count.get(slot) + 1);
+        }
+        self.sum.set(slot, self.sum.get(slot) + value);
+        if seen_before == 0 || value < self.min.get(slot) {
+            self.min.set(slot, value);
+        }
+        if value > self.max.get(slot) {
+            self.max.set(slot, value);
+        }
+    }
+
+    /// `min` needs a row implicitly missing a group member corrected to
+    /// `0` (see module docs); `presence_count` at `threshold == 0` is
+    /// trivially `group.len()` everywhere, computed directly rather than
+    /// trusting the scan (which never visits a cell that was never
+    /// present, so it would otherwise undercount).
+    fn finish(mut self, n: usize, group_len: u32, threshold: u32) -> io::Result<IntGroupStats> {
+        for slot in 0..n {
+            if self.present_count.get(slot) < group_len {
+                self.min.set(slot, 0);
+            }
+            if threshold == 0 {
+                self.presence_count.set(slot, group_len);
+            }
+        }
+        Ok(IntGroupStats {
+            presence_count: self.presence_count.freeze()?,
+            sum: self.sum.freeze()?,
+            min: self.min.freeze()?,
+            max: self.max.freeze()?,
+        })
+    }
+}
+
+/// See the module docs. `groups` may reference overlapping columns; each
+/// column is read once regardless of how many groups need it.
+pub fn batch_int_group_stats(
+    mat: &PersistentIntMatrix,
+    groups: &[ColGroup],
+    threshold: u32,
+) -> io::Result<Vec<IntGroupStats>> {
+    let n = mat.n();
+
+    let mut col_to_groups: HashMap<usize, Vec<usize>> = HashMap::new();
+    for (gi, g) in groups.iter().enumerate() {
+        for &c in &g.indices {
+            col_to_groups.entry(c).or_default().push(gi);
+        }
+    }
+
+    let mut accums: Vec<GroupAccum> = groups
+        .iter()
+        .map(|_| GroupAccum::new(n))
+        .collect::<io::Result<_>>()?;
+
+    match mat {
+        PersistentIntMatrix::Sparse(m) => {
+            for slot in 0..n {
+                m.for_each_cell_in_row(slot, |col, value| {
+                    if let Some(gs) = col_to_groups.get(&col) {
+                        for &gi in gs {
+                            accums[gi].touch(slot, value, threshold);
+                        }
+                    }
+                });
+            }
+        }
+        PersistentIntMatrix::Columnar(_) | PersistentIntMatrix::Packed(_) => {
+            for (&col, gs) in &col_to_groups {
+                let view = mat.col_view(col);
+                for (slot, value) in view.iter().enumerate() {
+                    if value == 0 {
+                        continue; // matches Sparse's own "never visit an absent cell"
+                    }
+                    for &gi in gs {
+                        accums[gi].touch(slot, value, threshold);
+                    }
+                }
+            }
+        }
+    }
+
+    accums
+        .into_iter()
+        .zip(groups)
+        .map(|(a, g)| a.finish(n, g.indices.len() as u32, threshold))
+        .collect()
+}
diff --git a/src/obicompactvec/src/lib.rs b/src/obicompactvec/src/lib.rs
index 3fbe47cd..fb1ad1ed 100644
--- a/src/obicompactvec/src/lib.rs
+++ b/src/obicompactvec/src/lib.rs
@@ -5,6 +5,7 @@ mod colgroup;
 mod eliasfano;
 mod fixedintvec;
 mod format;
+mod int_group_ops;
 mod intmatrix;
 mod layer_meta;
 mod matrix_builder;
@@ -27,6 +28,7 @@ pub use bitvec::{BitIter, PersistentBitVec, PersistentBitVecBuilder};
 pub use builder::PersistentCompactIntVecBuilder;
 pub use colgroup::{ColGroup, FilterMask, MatrixGroupOps, eval_filter_mask};
 pub use eliasfano::{EliasFano, EliasFanoBuilder};
+pub use int_group_ops::{IntGroupStats, batch_int_group_stats};
 pub use fixedintvec::{PersistentFixedIntVec, PersistentFixedIntVecBuilder, bit_width_for_range};
 pub use intmatrix::{
     PersistentCompactIntMatrixBuilder, PersistentIntMatrix, pack_compact_int_matrix,
diff --git a/src/obicompactvec/src/sparse_intmatrix.rs b/src/obicompactvec/src/sparse_intmatrix.rs
index 2d596f48..23675f4a 100644
--- a/src/obicompactvec/src/sparse_intmatrix.rs
+++ b/src/obicompactvec/src/sparse_intmatrix.rs
@@ -93,7 +93,7 @@ impl PersistentSparseCompactIntMatrix {
     /// ([`PersistentSparseBitMatrix::for_each_genome_in_row`]) with this
     /// row's private slice of whichever value stream it belongs to.
     #[inline]
-    fn for_each_cell_in_row(&self, slot: usize, mut f: impl FnMut(usize, u32)) {
+    pub(crate) fn for_each_cell_in_row(&self, slot: usize, mut f: impl FnMut(usize, u32)) {
         match self.support.row_rank(slot) {
             RowRank::Singleton(rank) => {
                 let v = self.singleton_values.get(rank);
diff --git a/src/obicompactvec/src/tests/int_group_ops.rs b/src/obicompactvec/src/tests/int_group_ops.rs
new file mode 100644
index 00000000..931af238
--- /dev/null
+++ b/src/obicompactvec/src/tests/int_group_ops.rs
@@ -0,0 +1,115 @@
+use tempfile::tempdir;
+
+use crate::{
+    ColGroup, PersistentCompactIntMatrixBuilder, PersistentIntMatrix,
+    PersistentSparseCompactIntMatrixBuilder, batch_int_group_stats,
+};
+
+fn make_columnar(cols: &[&[u32]]) -> (tempfile::TempDir, PersistentIntMatrix) {
+    let n = cols.first().map_or(0, |c| c.len());
+    let dir = tempdir().unwrap();
+    let mut b = PersistentCompactIntMatrixBuilder::new(n, &dir.path().join("counts")).unwrap();
+    for &col in cols {
+        let mut cb = b.add_col().unwrap();
+        for (slot, &v) in col.iter().enumerate() {
+            cb.set(slot, v);
+        }
+        cb.close().unwrap();
+    }
+    b.close().unwrap();
+    let m = PersistentIntMatrix::open(dir.path()).unwrap();
+    (dir, m)
+}
+
+/// Row-major input (`cols[c][row]`), sparse-encoded — zero values are
+/// never pushed, matching the format's own "nonzero cells only" contract.
+fn make_sparse(cols: &[&[u32]]) -> (tempfile::TempDir, PersistentIntMatrix) {
+    let n_cols = cols.len();
+    let n_rows = cols.first().map_or(0, |c| c.len());
+    let dir = tempdir().unwrap();
+    let sparse_dir = dir.path().join("counts");
+    let mut b = PersistentSparseCompactIntMatrixBuilder::new(n_rows, n_cols, &sparse_dir).unwrap();
+    let mut idx = Vec::new();
+    let mut vals = Vec::new();
+    for row in 0..n_rows {
+        idx.clear();
+        vals.clear();
+        for (c, col) in cols.iter().enumerate() {
+            if col[row] != 0 {
+                idx.push(c as u32);
+                vals.push(col[row]);
+            }
+        }
+        b.push_row(&idx, &vals);
+    }
+    let m = PersistentIntMatrix::Sparse(b.finish().unwrap());
+    (dir, m)
+}
+
+// col0=[3,0,5,2], col1=[0,4,5,0], col2=[1,4,0,2]
+// g1={0,1} (disjoint columns from g2's col2), g2={1,2} (shares col1 with g1)
+const COL0: [u32; 4] = [3, 0, 5, 2];
+const COL1: [u32; 4] = [0, 4, 5, 0];
+const COL2: [u32; 4] = [1, 4, 0, 2];
+
+#[test]
+fn batch_int_group_stats_matches_hand_computed_for_columnar_and_sparse() {
+    let cols: [&[u32]; 3] = [&COL0, &COL1, &COL2];
+    let (_dcol, columnar) = make_columnar(&cols);
+    let (_dsparse, sparse) = make_sparse(&cols);
+
+    let groups = [ColGroup::new("g1", vec![0, 1]), ColGroup::new("g2", vec![1, 2])];
+    let threshold = 2;
+
+    // g1={0,1}: presence_count, sum, min, max per row.
+    let g1_presence = [1u32, 1, 2, 1];
+    let g1_sum = [3u64, 4, 10, 2];
+    let g1_min = [0u32, 0, 5, 0]; // row2: both present (3>=thr? irrelevant to min) -> min(5,5)=5
+    let g1_max = [3u32, 4, 5, 2];
+
+    // g2={1,2}: col1 vs col2.
+    let g2_presence = [0u32, 2, 1, 1];
+    let g2_sum = [1u64, 8, 5, 2];
+    let g2_min = [0u32, 4, 0, 0];
+    let g2_max = [1u32, 4, 5, 2];
+
+    for (label, mat) in [("columnar", &columnar), ("sparse", &sparse)] {
+        let stats = batch_int_group_stats(mat, &groups, threshold).unwrap();
+        assert_eq!(stats.len(), 2, "{label}: wrong number of groups");
+
+        for row in 0..4 {
+            assert_eq!(stats[0].presence_count.get(row), g1_presence[row], "{label}: g1 presence row {row}");
+            assert_eq!(stats[0].sum.get(row) as u64, g1_sum[row], "{label}: g1 sum row {row}");
+            assert_eq!(stats[0].min.get(row), g1_min[row], "{label}: g1 min row {row}");
+            assert_eq!(stats[0].max.get(row), g1_max[row], "{label}: g1 max row {row}");
+
+            assert_eq!(stats[1].presence_count.get(row), g2_presence[row], "{label}: g2 presence row {row}");
+            assert_eq!(stats[1].sum.get(row) as u64, g2_sum[row], "{label}: g2 sum row {row}");
+            assert_eq!(stats[1].min.get(row), g2_min[row], "{label}: g2 min row {row}");
+            assert_eq!(stats[1].max.get(row), g2_max[row], "{label}: g2 max row {row}");
+        }
+    }
+}
+
+#[test]
+fn threshold_zero_makes_presence_count_the_group_size_everywhere() {
+    let cols: [&[u32]; 3] = [&COL0, &COL1, &COL2];
+    let (_dcol, columnar) = make_columnar(&cols);
+    let (_dsparse, sparse) = make_sparse(&cols);
+    let groups = [ColGroup::new("g1", vec![0, 1, 2])];
+
+    for (label, mat) in [("columnar", &columnar), ("sparse", &sparse)] {
+        let stats = batch_int_group_stats(mat, &groups, 0).unwrap();
+        for row in 0..4 {
+            assert_eq!(stats[0].presence_count.get(row), 3, "{label}: row {row}");
+        }
+    }
+}
+
+#[test]
+fn empty_group_list_returns_empty() {
+    let cols: [&[u32]; 1] = [&COL0];
+    let (_d, m) = make_columnar(&cols);
+    let stats = batch_int_group_stats(&m, &[], 1).unwrap();
+    assert!(stats.is_empty());
+}
diff --git a/src/obicompactvec/src/tests/mod.rs b/src/obicompactvec/src/tests/mod.rs
index 4f6fd9ac..ae7b1e5a 100644
--- a/src/obicompactvec/src/tests/mod.rs
+++ b/src/obicompactvec/src/tests/mod.rs
@@ -3,6 +3,7 @@ mod bitvec;
 mod colgroup;
 mod eliasfano;
 mod fixedintvec;
+mod int_group_ops;
 mod intmatrix;
 mod rankselect;
 mod sparse;
diff --git a/src/obikselect/src/select_layer.rs b/src/obikselect/src/select_layer.rs
index d2630ce0..9a668905 100644
--- a/src/obikselect/src/select_layer.rs
+++ b/src/obikselect/src/select_layer.rs
@@ -1,18 +1,24 @@
 //! Per-layer column projection/aggregation — the mechanics
-//! [`crate::select::Select`] runs once per partition-layer. Generic over
-//! source content (`Count`/`Presence`) via `obicompactvec::MatrixGroupOps`
-//! (object-safe: one `&dyn MatrixGroupOps` covers both matrix kinds, so
-//! this never branches on source content beyond the single `match` that
-//! opens it) and over destination content via [`DstBuilder`] — no
-//! duplicated per-op-per-content code path.
+//! [`crate::select::Select`] runs once per partition-layer. `Presence` and
+//! `Count` sources each get their own batched, row-major-friendly pass
+//! (`obicompactvec::batch_presence_counts`/`batch_int_group_stats`)
+//! computing every output spec's statistics together in one shared scan,
+//! rather than the generic `MatrixGroupOps` per-group dispatch this used to
+//! go through (still used by `obikfilter`, whose ad hoc `FilterMask`
+//! column-index lists don't fit the same "one shared batch" shape) —
+//! `sum`/`min`/`max` mean genuinely different things for the two content
+//! kinds (a bit matrix's are cheap derivations of one presence count, a
+//! count matrix's are real per-value reductions), so each gets its own
+//! `agg_result_from_*` derivation, not a shared `AggOp` dispatch.
 
 use std::fs;
 use std::io;
 use std::path::Path;
 
 use obicompactvec::{
-    ColGroup, MatrixBuilder, MatrixGroupOps, PersistentBitMatrix, PersistentIntMatrix, TempBitVec,
-    TempBitVecBuilder, TempCompactIntVec, TempCompactIntVecBuilder, batch_presence_counts,
+    ColGroup, IntGroupStats, MatrixBuilder, PersistentBitMatrix, PersistentIntMatrix, TempBitVec,
+    TempBitVecBuilder, TempCompactIntVec, TempCompactIntVecBuilder, batch_int_group_stats,
+    batch_presence_counts,
 };
 use obikindex::layer::{KmerLayer, LayerContent};
 use obikindex::{KmerIndex, OKIError, OKIResult};
@@ -77,22 +83,6 @@ enum AggResult {
     Int(TempCompactIntVec),
 }
 
-fn compute_group(
-    mat: &dyn MatrixGroupOps,
-    spec: &OutputCol,
-    threshold: u32,
-) -> io::Result<AggResult> {
-    let g = ColGroup::new(spec.label.clone(), spec.indices.clone());
-    Ok(match spec.op {
-        AggOp::Any => AggResult::Bit(mat.partial_group_any(&g, threshold)?),
-        AggOp::All => AggResult::Bit(mat.partial_group_all(&g, threshold)?),
-        AggOp::None => AggResult::Bit(mat.partial_group_none(&g, threshold)?),
-        AggOp::Sum => AggResult::Int(mat.partial_group_sum(&g)?),
-        AggOp::Min => AggResult::Int(mat.partial_group_min(&g)?),
-        AggOp::Max => AggResult::Int(mat.partial_group_max(&g)?),
-    })
-}
-
 /// A bit-matrix source's every `AggOp` is a cheap derivation of one shared
 /// presence count (see [`obicompactvec::batch_presence_counts`]'s own
 /// docs for why: `sum` = the count itself, `any`/`max` = `count ≥ 1`,
@@ -131,6 +121,35 @@ fn agg_result_from_count(op: AggOp, group_len: usize, count: TempCompactIntVec)
     })
 }
 
+/// A count-matrix source's `AggOp` from [`IntGroupStats`] — unlike the
+/// bit-matrix case, `sum`/`min`/`max` are genuine per-value reductions
+/// already computed by `batch_int_group_stats`, not further derivations;
+/// only `any`/`all`/`none` still need a threshold-count comparison here.
+fn agg_result_from_int_stats(op: AggOp, group_len: usize, stats: IntGroupStats) -> io::Result<AggResult> {
+    let n = stats.presence_count.len();
+    let group_len = group_len as u32;
+    Ok(match op {
+        AggOp::Sum => AggResult::Int(stats.sum),
+        AggOp::Min => AggResult::Int(stats.min),
+        AggOp::Max => AggResult::Int(stats.max),
+        AggOp::Any => {
+            let mut b = TempBitVecBuilder::new(n)?;
+            b.or_where(stats.presence_count.view(), |v| v >= 1);
+            AggResult::Bit(b.freeze()?)
+        }
+        AggOp::All => {
+            let mut b = TempBitVecBuilder::new(n)?;
+            b.or_where(stats.presence_count.view(), |v| v == group_len);
+            AggResult::Bit(b.freeze()?)
+        }
+        AggOp::None => {
+            let mut b = TempBitVecBuilder::new(n)?;
+            b.or_where(stats.presence_count.view(), |v| v == 0);
+            AggResult::Bit(b.freeze()?)
+        }
+    })
+}
+
 // ── AggResult → MatrixBuilder ─────────────────────────────────────────────────
 
 /// Add one already-aggregated column to `mb` — the only piece `MatrixBuilder`
@@ -233,8 +252,14 @@ pub(crate) fn select_partition(
             // unchanged, one `col_view`-driven pass per spec.
             LayerContent::Count => {
                 let mat = PersistentIntMatrix::open(src_layer.dir()).map_err(OKIError::Io)?;
-                for spec in specs {
-                    let r = compute_group(&mat, spec, threshold).map_err(OKIError::Io)?;
+                let groups: Vec<ColGroup> = specs
+                    .iter()
+                    .map(|s| ColGroup::new(s.label.clone(), s.indices.clone()))
+                    .collect();
+                let stats = batch_int_group_stats(&mat, &groups, threshold).map_err(OKIError::Io)?;
+                for (spec, stat) in specs.iter().zip(stats) {
+                    let r = agg_result_from_int_stats(spec.op, spec.indices.len(), stat)
+                        .map_err(OKIError::Io)?;
                     add_result(&mut builder, r).map_err(OKIError::Io)?;
                 }
             }