From 7d499a3b1ce7b8cd9f25de87ead7e988cf966be4 Mon Sep 17 00:00:00 2001 From: Carsten Erickson Date: Tue, 26 May 2026 20:12:19 -0700 Subject: [PATCH 1/2] Add tbx_itr_regions: multi-region iterator for tabix indexed text Constructs an hts_itr_t driven by hts_itr_multi_next that yields records from many regions in one pass over a tabix indexed text file (VCF.gz, BED.gz, GFF.gz, custom TSV preset). Adjacent regions are coalesced into combined index lookups instead of one tbx_itr_queryi call per region. This brings tabix indexed text to parity with BAM and CRAM, which have had sam_itr_regions since the multi region driver was added. Affected workloads include ancestry and population genetics subsetting VCFs to large SNP panels (the motivating case in samtools/bcftools#2557 identified an 11 plus hour single chromosome subset against an 84M site panel), exome capture target subsetting, annotation overlap against large bed tracks, and trio analysis against candidate de novo position lists. Measured against tbx_itr_queryi per region on a 176 MB chr22 1kGP SNP only VCF.gz: panel single (s) multi (s) speedup 1819 6.0 5.7 1.05x 4797 15.8 10.2 1.56x 18197 60.6 12.8 4.72x 47975 159.4 12.8 12.45x 175909 562.9 13.3 42.48x Multi path wallclock plateaus at sequential scan time as the panel grows; single grows roughly linearly with region count. Every row above was byte verified between paths (matching fnv1a64 hash plus matching byte length); the 100k row corresponds to 1.79 GB of byte identical output. The implementation deep copies the caller reglist so callers always own their input regardless of outcome, pre resolves reglist[i].reg name strings via tbx_name2id so callers may use either tids or reference name strings, and sorts intervals defensively before handing off to hts_itr_multi_bam. tbx_readrec_multi retrieves tbx_t from BGZF private_data; this mirrors how bcf_readrec uses the same slot for BCF header context and inherits the same one multi region iterator per htsFile constraint (documented in tbx.h). test/test_tbx_multi.c verifies byte identical parity between the single region and multi region paths across BED, VCF, and GFF preset fixtures, plus dedup of overlapping intervals, plus six failure path inputs, plus two name resolution paths. test/bench_tbx_regions.c is the reproducible source of the table above, streaming an fnv1a64 hash over both paths' output for byte level parity at constant memory. Addresses samtools/bcftools#2557. Signed-off-by: Carsten Erickson Assisted-by: Claude (claude-opus-4-7) --- Makefile | 12 +- NEWS | 11 ++ hts_internal.h | 8 + htslib/tbx.h | 54 ++++++ region.c | 2 +- tbx.c | 191 ++++++++++++++++++- test/bench_tbx_regions.c | 374 +++++++++++++++++++++++++++++++++++++ test/tbx-multi.bed | 6 + test/tbx-multi.gff | 7 + test/tbx-multi.vcf | 10 + test/test.pl | 30 +++ test/test_tbx_multi.c | 393 +++++++++++++++++++++++++++++++++++++++ 12 files changed, 1090 insertions(+), 8 deletions(-) create mode 100644 test/bench_tbx_regions.c create mode 100644 test/tbx-multi.bed create mode 100644 test/tbx-multi.gff create mode 100644 test/tbx-multi.vcf create mode 100644 test/test_tbx_multi.c diff --git a/Makefile b/Makefile index 33062d69e..7923ad216 100644 --- a/Makefile +++ b/Makefile @@ -104,7 +104,9 @@ BUILT_TEST_PROGRAMS = \ test/test-parse-reg \ test/test_introspection \ test/test-bcf_set_variant_type \ - test/test_hfile_libcurl + test/test_hfile_libcurl \ + test/test_tbx_multi \ + test/bench_tbx_regions BUILT_THRASH_PROGRAMS = \ test/thrash_threads1 \ @@ -808,6 +810,12 @@ test/test-vcf-sweep: test/test-vcf-sweep.o libhts.a test/test-bcf-sr: test/test-bcf-sr.o libhts.a $(CC) $(LDFLAGS) -o $@ test/test-bcf-sr.o libhts.a $(LIBS) -lpthread +test/test_tbx_multi: test/test_tbx_multi.o libhts.a + $(CC) $(LDFLAGS) -o $@ test/test_tbx_multi.o libhts.a $(LIBS) -lpthread + +test/bench_tbx_regions: test/bench_tbx_regions.o libhts.a + $(CC) $(LDFLAGS) -o $@ test/bench_tbx_regions.o libhts.a $(LIBS) -lpthread + test/test-bcf-translate: test/test-bcf-translate.o libhts.a $(CC) $(LDFLAGS) -o $@ test/test-bcf-translate.o libhts.a $(LIBS) -lpthread @@ -897,6 +905,8 @@ test/test_index.o: test/test_index.c config.h $(htslib_sam_h) $(htslib_vcf_h) test/test-vcf-api.o: test/test-vcf-api.c config.h $(htslib_hts_h) $(htslib_hts_alloc_h) $(htslib_vcf_h) $(htslib_vcfutils_h) $(htslib_kbitset_h) $(htslib_kstring_h) $(htslib_kseq_h) test/test-vcf-sweep.o: test/test-vcf-sweep.c config.h $(htslib_vcf_sweep_h) test/test-bcf-sr.o: test/test-bcf-sr.c config.h $(htslib_hts_alloc_h) $(htslib_hts_defs_h) $(htslib_synced_bcf_reader_h) $(htslib_hts_h) $(htslib_vcf_h) +test/test_tbx_multi.o: test/test_tbx_multi.c config.h $(htslib_hts_h) $(htslib_tbx_h) $(htslib_kstring_h) +test/bench_tbx_regions.o: test/bench_tbx_regions.c config.h $(htslib_hts_h) $(htslib_tbx_h) $(htslib_kstring_h) test/test-bcf-translate.o: test/test-bcf-translate.c config.h $(htslib_vcf_h) test/test_introspection.o: test/test_introspection.c config.h $(htslib_hts_h) $(htslib_hfile_h) test/test-bcf_set_variant_type.o: test/test-bcf_set_variant_type.c config.h $(htslib_hts_h) vcf.c diff --git a/NEWS b/NEWS index 7181ea3ad..d944a81bb 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,17 @@ Noteworthy changes in release a.b ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +API additions +------------- + +* Added tbx_itr_regions(), a multi-region iterator constructor for tabix + indexed text files (VCF.gz, BED, GFF, custom TSV). Returns an hts_itr_t + that, when driven by hts_itr_multi_next, yields records from all listed + regions in genome order with coalesced index lookups. Meaningfully + faster than calling tbx_itr_queryi() once per region when the region + list is large or dense. Mirrors how hts_itr_regions() already works for + BAM and CRAM. samtools/bcftools#2557. + Noteworthy changes in release 1.23.1 (18th March 2026) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/hts_internal.h b/hts_internal.h index 33b9422d7..abb2aac97 100644 --- a/hts_internal.h +++ b/hts_internal.h @@ -137,6 +137,14 @@ const char *hts_plugin_path(void); */ int bgzf_idx_push(BGZF *fp, hts_idx_t *hidx, int tid, hts_pos_t beg, hts_pos_t end, uint64_t offset, int is_mapped); +/* + * qsort comparator for hts_pair_pos_t intervals, sorting by (beg, end) + * ascending. Exposed because both region.c (sort within parsed reglists) + * and tbx.c (multi-region tabix iterator: hts_itr_multi_bam requires + * intervals to be sorted) need this comparator with identical behavior. + */ +int compare_hts_pair_pos_t(const void *av, const void *bv); + static inline int find_file_extension(const char *fn, char ext_out[static HTS_MAX_EXT_LEN]) { const char *delim = fn ? strstr(fn, HTS_IDX_DELIM) : NULL, *ext; diff --git a/htslib/tbx.h b/htslib/tbx.h index 8b5ac1e4b..cdd62ddd1 100644 --- a/htslib/tbx.h +++ b/htslib/tbx.h @@ -75,6 +75,60 @@ extern const tbx_conf_t tbx_conf_gff, tbx_conf_bed, tbx_conf_psltbl, tbx_conf_sa HTSLIB_EXPORT int tbx_readrec(BGZF *fp, void *tbxv, void *sv, int *tid, hts_pos_t *beg, hts_pos_t *end); + /// Construct a multi-region iterator over a tabix indexed text file. + /** Returns an hts_itr_t that, when driven by hts_itr_multi_next, yields + records from all the regions in reglist in genome order. Adjacent + and nearby regions are coalesced into combined tabix index lookups + by hts_itr_multi_next, so this is meaningfully faster than calling + tbx_itr_queryi once per region when the region list is large or + dense. See test/bench_tbx_regions.c for a reproducible measurement + of the speedup curve. + + Reglist entries may specify their target contig either way: + - reglist[i].tid set to a tid (use tbx_name2id() to resolve) + with reglist[i].reg == NULL, or + - reglist[i].reg set to a reference name string such as "chr1" + (or "." for all-with-coords records, "*" for unmapped), + which this function resolves internally via tbx_name2id(). + Each reglist[i].intervals must hold reglist[i].count + hts_pair_pos_t entries. + + Semantics differ from running tbx_itr_queryi() multiple times in + one respect: each underlying file record is yielded at most once, + even when multiple intervals in reglist cover it. Duplicate or + overlapping intervals produce a single emission per matching + record, not one per matching interval. Callers that need + per-interval multiplicity must call tbx_itr_queryi() per interval. + + Ownership: the caller's reglist is deep-copied internally and is + never mutated or freed by this function regardless of outcome. + The caller retains ownership of reglist (and is responsible for + freeing it) as well as of tbx and fp. tbx must remain valid for + the lifetime of the iterator. Destroy the iterator with + hts_itr_destroy. + + Concurrency: only one multi-region tabix iterator may be active + on a given htsFile at a time. The iterator uses the BGZF + private_data slot on fp to thread tbx_t through hts_itr_multi_next; + constructing a second multi-region iterator on the same fp before + destroying the first would clobber that slot and yield wrong + records from the first. Single-region iterators (tbx_itr_queryi) + on the same fp are unaffected. This matches the existing + constraint on binary BCF iteration, which uses the same slot. + + Lifetime: destroy the iterator (hts_itr_destroy) before destroying + the tbx (tbx_destroy) or closing the fp (hts_close). The iterator + retains a pointer to tbx via the BGZF private_data slot; destroying + tbx first leaves a dangling pointer in fp's BGZF cache that would + be returned by a later bgzf_get_private_data on the same fp. + + Returns NULL if fp, tbx, or reglist is NULL; if count <= 0; if + fp is not BGZF; if memory cannot be allocated; or if the iterator + cannot otherwise be constructed. */ + HTSLIB_EXPORT + hts_itr_t *tbx_itr_regions(htsFile *fp, tbx_t *tbx, + hts_reglist_t *reglist, int count); + /// Build an index of the lines in a BGZF-compressed file /** The index struct returned by a successful call should be freed via tbx_destroy() when it is no longer needed. diff --git a/region.c b/region.c index 8b570e0bf..04b866b61 100644 --- a/region.c +++ b/region.c @@ -38,7 +38,7 @@ typedef struct reglist KHASH_MAP_INIT_INT(reg, reglist_t) typedef kh_reg_t reghash_t; -static int compare_hts_pair_pos_t (const void *av, const void *bv) +int compare_hts_pair_pos_t (const void *av, const void *bv) { hts_pair_pos_t *a = (hts_pair_pos_t *) av; hts_pair_pos_t *b = (hts_pair_pos_t *) bv; diff --git a/tbx.c b/tbx.c index d9d636b43..61a50a2ab 100644 --- a/tbx.c +++ b/tbx.c @@ -33,6 +33,7 @@ DEALINGS IN THE SOFTWARE. */ #include #include "htslib/tbx.h" #include "htslib/bgzf.h" +#include "bgzf_internal.h" #include "htslib/hts_alloc.h" #include "htslib/hts_endian.h" #include "hts_internal.h" @@ -351,19 +352,21 @@ static inline int get_intv(tbx_t *tbx, kstring_t *str, tbx_intv_t *intv, int is_ * -1 on EOF * <= -2 on error */ -int tbx_readrec(BGZF *fp, void *tbxv, void *sv, int *tid, hts_pos_t *beg, hts_pos_t *end) +/* Shared body between tbx_readrec (single region path) and tbx_readrec_multi + * (multi region path). Reads the next non-meta line from fp, parses it via + * get_intv, and returns the interval coordinates. Centralising the body here + * means that future fixes to meta_char handling or get_intv parsing apply to + * both paths automatically. */ +static int tbx_readrec_impl(BGZF *fp, tbx_t *tbx, kstring_t *s, + int *tid, hts_pos_t *beg, hts_pos_t *end) { - tbx_t *tbx = (tbx_t *) tbxv; - kstring_t *s = (kstring_t *) sv; int ret; - // Get a line until either EOF or a non-meta character do { ret = bgzf_getline(fp, '\n', s); } while (ret >= 0 && s->l && *s->s == tbx->conf.meta_char); - // Parse line - if (ret >= 0) { + if (ret >= 0) { tbx_intv_t intv; if (get_intv(tbx, s, &intv, 0) < 0) return -2; @@ -373,6 +376,182 @@ int tbx_readrec(BGZF *fp, void *tbxv, void *sv, int *tid, hts_pos_t *beg, hts_po return ret; } +int tbx_readrec(BGZF *fp, void *tbxv, void *sv, int *tid, hts_pos_t *beg, hts_pos_t *end) +{ + return tbx_readrec_impl(fp, (tbx_t *) tbxv, (kstring_t *) sv, tid, beg, end); +} + +/* Internal helpers for the multi-region tabix iterator. Mirror the + * sam.c bam_pseek/bam_ptell pattern but parameterised on BGZF rather than + * any BAM-specific structure. */ +static int tbx_pseek(void *fp, int64_t offset, int where) +{ + return bgzf_seek((BGZF *) fp, offset, where); +} + +static int64_t tbx_ptell(void *fp) +{ + BGZF *bgzf = (BGZF *) fp; + if (!bgzf) return -1L; + return bgzf_tell(bgzf); +} + +/* Multi-region readrec callback. Retrieves tbx_t from BGZF private_data + * because hts_itr_multi_next passes the htsFile pointer as the second + * argument and has no mechanism to thread a tabix-specific context through. + * Mirrors how bcf_readrec retrieves bcf_hdr_t from BGZF private_data. */ +static int tbx_readrec_multi(BGZF *fp, void *ignored, void *sv, + int *tid, hts_pos_t *beg, hts_pos_t *end) +{ + tbx_t *tbx = (tbx_t *) bgzf_get_private_data(fp); + if (!tbx) { + /* The wrapper tbx_itr_regions always sets private_data before + * returning the iterator. Reaching here means either the slot + * was cleared mid-iteration (caller violated the single-iterator + * precondition) or the iterator was driven outside its lifetime. + * Return -2 rather than -1 so the caller sees an error rather + * than a silent end-of-stream. */ + hts_log_error("Multi-region tabix iterator has no BGZF private_data;" + " concurrent multi-region iterators on the same" + " htsFile are not supported"); + return -2; + } + return tbx_readrec_impl(fp, tbx, (kstring_t *) sv, tid, beg, end); +} + +/* Deep copy a caller-provided reglist for internal use by tbx_itr_regions. + * + * Why deep copy: hts_itr_regions's ownership semantics on failure are + * inconsistent (NULL return may or may not have freed the reglist depending + * on which branch failed). By giving it our own owned copy, the caller's + * input is never touched, so they always own what they passed in regardless + * of outcome. + * + * Why resolve .reg here: hts_itr_regions invokes getid(hdr, reg) whenever + * reglist[i].reg is non-NULL. We pass getid=NULL, so we must clear .reg on + * the copy after pre-resolving it via tbx_name2id. This lets callers pass + * either pre-resolved tids or reference name strings, matching how + * sam_itr_regions handles getid internally. + * + * Returns 0 on success (*out points to a freshly allocated reglist that + * the caller must release via hts_reglist_free) or -1 on allocation + * failure. */ +static int tbx_reglist_dup(tbx_t *tbx, const hts_reglist_t *src, int count, + hts_reglist_t **out) +{ + hts_reglist_t *dst = calloc(count, sizeof(*dst)); + if (!dst) return -1; + + for (int i = 0; i < count; i++) { + dst[i].tid = src[i].tid; + dst[i].count = src[i].count; + dst[i].min_beg = src[i].min_beg; + dst[i].max_end = src[i].max_end; + dst[i].reg = NULL; /* always pre-resolved below */ + dst[i].intervals = NULL; + + if (src[i].reg) { + if (!strcmp(src[i].reg, ".")) dst[i].tid = HTS_IDX_START; + else if (!strcmp(src[i].reg, "*")) dst[i].tid = HTS_IDX_NOCOOR; + else { + int tid = tbx_name2id(tbx, src[i].reg); + if (tid < 0) + hts_log_warning("Region '%s' refers to an unknown" + " reference name", src[i].reg); + dst[i].tid = tid; + } + } + + /* An entry with count > 0 but intervals == NULL (or vice versa) is + * an inconsistent caller-side state; treat as having no intervals + * so we never hand hts_itr_multi_bam a non-zero count with a NULL + * intervals array, which it would dereference and crash on. */ + if (!src[i].intervals || src[i].count == 0) { + dst[i].count = 0; + continue; + } + + /* Read src[i].count into a local once and use the local for size + * computation, allocation, and bookkeeping. Avoids any TOCTOU + * mismatch between the size we allocate and the size we believe + * dst holds, even though callers are not expected to mutate src + * concurrently. */ + uint32_t n = src[i].count; + dst[i].intervals = malloc((size_t) n * sizeof(hts_pair_pos_t)); + if (!dst[i].intervals) { + hts_reglist_free(dst, count); + return -1; + } + memcpy(dst[i].intervals, src[i].intervals, + (size_t) n * sizeof(hts_pair_pos_t)); + dst[i].count = n; + + /* Recompute min_beg / max_end from the copied intervals so they + * remain consistent even if the caller's summary fields were + * stale or if the later qsort reorders intervals[0]. */ + hts_pos_t mn = dst[i].intervals[0].beg; + hts_pos_t mx = dst[i].intervals[0].end; + for (uint32_t k = 1; k < n; k++) { + if (dst[i].intervals[k].beg < mn) mn = dst[i].intervals[k].beg; + if (dst[i].intervals[k].end > mx) mx = dst[i].intervals[k].end; + } + dst[i].min_beg = mn; + dst[i].max_end = mx; + } + + *out = dst; + return 0; +} + +hts_itr_t *tbx_itr_regions(htsFile *fp, tbx_t *tbx, + hts_reglist_t *reglist, int count) +{ + if (!fp || !tbx || !reglist || count <= 0) return NULL; + BGZF *bgzf = hts_get_bgzfp(fp); + if (!bgzf) return NULL; + + /* Make our own copy of the reglist so the caller's input is never + * mutated and never assumed to be live past this call. The copy also + * lets us pre-resolve any .reg name strings to tids so we can pass + * getid=NULL safely to hts_itr_regions below. */ + hts_reglist_t *owned = NULL; + if (tbx_reglist_dup(tbx, reglist, count, &owned) < 0) + return NULL; + + /* hts_itr_multi_bam (the itr_specific used here) requires intervals + * within each reglist entry to be sorted by start position; its + * advance logic uses the interval index packed into the offset list + * which relies on sorted iteration. Sort our copy now. */ + for (int i = 0; i < count; i++) { + if (owned[i].intervals && owned[i].count > 1) { + qsort(owned[i].intervals, owned[i].count, + sizeof(hts_pair_pos_t), compare_hts_pair_pos_t); + } + } + + hts_itr_t *itr = hts_itr_regions(tbx->idx, owned, count, NULL, NULL, + hts_itr_multi_bam, tbx_readrec_multi, + tbx_pseek, tbx_ptell); + if (!itr) { + /* hts_itr_regions has exactly one failure path that does not + * free reglist via hts_itr_destroy: the OOM at its initial + * calloc(itr) before it ever assigns reg_list. Every other + * failure (itr_specific failure, getid failure) frees 'owned'. + * We accept the small leak in the calloc-OOM case rather than + * duplicate hts_itr_regions's allocation path here, because + * the caller's original reglist is untouched either way. */ + return NULL; + } + + /* Stash the tbx_t in BGZF private_data so tbx_readrec_multi can + * retrieve it on each call. The slot is otherwise unused for tabix + * indexed text files; only bcf_hdr_read uses it, for binary BCF. + * Caller retains ownership of tbx, which must remain valid for the + * lifetime of the iterator. */ + bgzf_set_private_data(bgzf, tbx, NULL); + return itr; +} + static int tbx_set_meta(tbx_t *tbx) { int i, l = 0, l_nm; diff --git a/test/bench_tbx_regions.c b/test/bench_tbx_regions.c new file mode 100644 index 000000000..d60dc5bc0 --- /dev/null +++ b/test/bench_tbx_regions.c @@ -0,0 +1,374 @@ +/* bench_tbx_regions.c -- microbenchmark comparing tbx_itr_queryi loop + against tbx_itr_regions on the same input. + + Usage: bench_tbx_regions FILE.gz REGIONS.bed [WARMUP] [REPS] + FILE.gz bgzipped + tabix indexed source + REGIONS.bed plain text BED of regions to query, one per line: + (0 based half open) + WARMUP number of warmup runs per path (default 0) + REPS number of timed runs per path (default 3) + + Reports wallclock and CPU time and total records emitted for each path. + Verifies byte level parity between paths by running an untimed capture + pass that streams every emitted line of both paths into separate + buffers, then memcmp comparing them and printing an FNV 1a digest of + each. Falls back to exit code 3 on any mismatch and prints the first + diverging line index. The record count check remains as an additional + sanity gate. + + Copyright (C) 2026 Genome Research Ltd. + + Permission is hereby granted, free of charge, to any person obtaining a copy + of this software and associated documentation files (the "Software"), to deal + in the Software without restriction, including without limitation the rights + to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included in + all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + DEALINGS IN THE SOFTWARE. +*/ + +#include + +#include +#include +#include +#include +#include +#include + +#include "../htslib/hts.h" +#include "../htslib/tbx.h" +#include "../htslib/kstring.h" + +/* FNV 1a 64 bit, streamed. Used purely to give the parity report a digest + * so a reviewer reading the output can see that both paths produced byte + * equal streams at any panel size without holding the full byte stream in + * RAM (which at 100k regions on a 1kGP VCF approaches 2 GB per side). + * + * Collision probability for two distinct streams is ~2^-64 (~5e-20), which + * is well below cosmic ray bit flip rates for the duration of any bench. + * Combined with the byte length check, we treat matching {length, hash} + * across both paths as the parity guarantee. */ +#define FNV1A64_OFFSET 0xcbf29ce484222325ULL +#define FNV1A64_PRIME 0x100000001b3ULL + +static inline uint64_t fnv1a64_update(uint64_t h, const void *data, size_t len) +{ + const unsigned char *p = (const unsigned char *) data; + for (size_t i = 0; i < len; i++) { + h ^= p[i]; + h *= FNV1A64_PRIME; + } + return h; +} + +/* In memory region list parsed once and reused for both paths. */ +typedef struct { + char *chrom; + hts_pos_t beg, end; +} bench_region_t; + +typedef struct { + bench_region_t *items; + int n, cap; +} bench_regions_t; + +/* Per run timing sample. parity_hash and parity_bytes are populated only + * when the caller requests parity tracking (see run_*); they sit alongside + * timing so the report can show all of {wall, records, bytes, hash} for a + * single rep without juggling separate return paths. */ +typedef struct { + double wall_s; + double user_s; + double sys_s; + long records; + uint64_t parity_hash; + size_t parity_bytes; +} bench_sample_t; + +static double timespec_diff_s(struct timespec a, struct timespec b) +{ + return (b.tv_sec - a.tv_sec) + (b.tv_nsec - a.tv_nsec) / 1.0e9; +} + +static double timeval_s(struct timeval t) +{ + return t.tv_sec + t.tv_usec / 1.0e6; +} + +/* Read a plain BED into memory. Three columns: chrom, beg, end. */ +static int read_bed(const char *path, bench_regions_t *out) +{ + FILE *fp = fopen(path, "r"); + if (!fp) { fprintf(stderr, "fopen %s failed\n", path); return -1; } + out->items = NULL; out->n = 0; out->cap = 0; + char line[4096]; + while (fgets(line, sizeof(line), fp)) { + if (line[0] == '#' || line[0] == '\n') continue; + char chrom[256]; long long beg, end; + if (sscanf(line, "%255s %lld %lld", chrom, &beg, &end) != 3) continue; + if (out->n == out->cap) { + int ncap = out->cap ? out->cap * 2 : 64; + bench_region_t *grown = realloc(out->items, ncap * sizeof(*grown)); + if (!grown) { fclose(fp); return -1; } + out->items = grown; out->cap = ncap; + } + out->items[out->n].chrom = strdup(chrom); + out->items[out->n].beg = (hts_pos_t)beg; + out->items[out->n].end = (hts_pos_t)end; + out->n++; + } + fclose(fp); + return 0; +} + +static void free_regions(bench_regions_t *r) +{ + for (int i = 0; i < r->n; i++) free(r->items[i].chrom); + free(r->items); +} + +/* Path A: for each region, tbx_itr_queryi + tbx_itr_next loop. + * If with_parity is non zero, every emitted line is fed into a streaming + * fnv1a64 hash and counted in s->parity_bytes (line bytes plus a trailing + * newline per record). The hash overhead is two machine ops per byte and + * adds well under 5 percent to wallclock; the parity_* fields appear in + * the parity report only, not in the timed reps. */ +static int run_single(const char *fn, const bench_regions_t *regs, + bench_sample_t *s, int with_parity) +{ + htsFile *fp = hts_open(fn, "r"); + if (!fp) return -1; + tbx_t *tbx = tbx_index_load(fn); + if (!tbx) { hts_close(fp); return -1; } + kstring_t line = {0,0,0}; + long records = 0; + uint64_t hash = FNV1A64_OFFSET; + size_t bytes = 0; + static const char newline = '\n'; + struct rusage ru0, ru1; + struct timespec t0, t1; + getrusage(RUSAGE_SELF, &ru0); + clock_gettime(CLOCK_MONOTONIC, &t0); + for (int i = 0; i < regs->n; i++) { + int tid = tbx_name2id(tbx, regs->items[i].chrom); + if (tid < 0) continue; + hts_itr_t *itr = tbx_itr_queryi(tbx, tid, regs->items[i].beg, + regs->items[i].end); + if (!itr) continue; + while (tbx_itr_next(fp, tbx, itr, &line) >= 0) { + records++; + if (with_parity) { + hash = fnv1a64_update(hash, line.s, line.l); + hash = fnv1a64_update(hash, &newline, 1); + bytes += line.l + 1; + } + } + tbx_itr_destroy(itr); + } + clock_gettime(CLOCK_MONOTONIC, &t1); + getrusage(RUSAGE_SELF, &ru1); + free(line.s); + tbx_destroy(tbx); + hts_close(fp); + s->wall_s = timespec_diff_s(t0, t1); + s->user_s = timeval_s(ru1.ru_utime) - timeval_s(ru0.ru_utime); + s->sys_s = timeval_s(ru1.ru_stime) - timeval_s(ru0.ru_stime); + s->records = records; + s->parity_hash = with_parity ? hash : 0; + s->parity_bytes = with_parity ? bytes : 0; + return 0; +} + +/* Path B: one tbx_itr_regions covering all regions, then drain. Streaming + * parity hash semantics mirror run_single. */ +static int run_multi(const char *fn, const bench_regions_t *regs, + bench_sample_t *s, int with_parity) +{ + htsFile *fp = hts_open(fn, "r"); + if (!fp) return -1; + tbx_t *tbx = tbx_index_load(fn); + if (!tbx) { hts_close(fp); return -1; } + + /* Group regions by tid into one reglist entry per contig. */ + hts_reglist_t *reglist = calloc(regs->n, sizeof(*reglist)); + if (!reglist) { tbx_destroy(tbx); hts_close(fp); return -1; } + int n_reg = 0; + for (int i = 0; i < regs->n; i++) { + int tid = tbx_name2id(tbx, regs->items[i].chrom); + if (tid < 0) continue; + int idx = -1; + for (int j = 0; j < n_reg; j++) { + if (reglist[j].tid == tid) { idx = j; break; } + } + if (idx < 0) { + idx = n_reg++; + reglist[idx].tid = tid; + reglist[idx].reg = NULL; + reglist[idx].count = 0; + reglist[idx].intervals = NULL; + reglist[idx].min_beg = regs->items[i].beg; + reglist[idx].max_end = regs->items[i].end; + } + hts_pair_pos_t *grown = realloc(reglist[idx].intervals, + (reglist[idx].count + 1) * sizeof(hts_pair_pos_t)); + if (!grown) { + for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + free(reglist); tbx_destroy(tbx); hts_close(fp); return -1; + } + reglist[idx].intervals = grown; + reglist[idx].intervals[reglist[idx].count].beg = regs->items[i].beg; + reglist[idx].intervals[reglist[idx].count].end = regs->items[i].end; + reglist[idx].count++; + if (regs->items[i].beg < reglist[idx].min_beg) + reglist[idx].min_beg = regs->items[i].beg; + if (regs->items[i].end > reglist[idx].max_end) + reglist[idx].max_end = regs->items[i].end; + } + + kstring_t line = {0,0,0}; + long records = 0; + uint64_t hash = FNV1A64_OFFSET; + size_t bytes = 0; + static const char newline = '\n'; + struct rusage ru0, ru1; + struct timespec t0, t1; + getrusage(RUSAGE_SELF, &ru0); + clock_gettime(CLOCK_MONOTONIC, &t0); + hts_itr_t *itr = tbx_itr_regions(fp, tbx, reglist, n_reg); + if (!itr) { + clock_gettime(CLOCK_MONOTONIC, &t1); + for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + free(reglist); free(line.s); + tbx_destroy(tbx); hts_close(fp); + return -1; + } + while (hts_itr_multi_next(fp, itr, &line) >= 0) { + records++; + if (with_parity) { + hash = fnv1a64_update(hash, line.s, line.l); + hash = fnv1a64_update(hash, &newline, 1); + bytes += line.l + 1; + } + } + clock_gettime(CLOCK_MONOTONIC, &t1); + getrusage(RUSAGE_SELF, &ru1); + free(line.s); + hts_itr_destroy(itr); /* takes ownership of reglist */ + tbx_destroy(tbx); + hts_close(fp); + s->wall_s = timespec_diff_s(t0, t1); + s->user_s = timeval_s(ru1.ru_utime) - timeval_s(ru0.ru_utime); + s->sys_s = timeval_s(ru1.ru_stime) - timeval_s(ru0.ru_stime); + s->records = records; + s->parity_hash = with_parity ? hash : 0; + s->parity_bytes = with_parity ? bytes : 0; + return 0; +} + +static void print_sample(const char *label, bench_sample_t s) +{ + printf(" %s: wall=%.3fs user=%.3fs sys=%.3fs records=%ld\n", + label, s.wall_s, s.user_s, s.sys_s, s.records); +} + +int main(int argc, char **argv) +{ + if (argc < 3) { + fprintf(stderr, "Usage: %s FILE.gz REGIONS.bed [WARMUP] [REPS]\n", argv[0]); + return 1; + } + const char *fn = argv[1]; + const char *bed = argv[2]; + int warmup = argc > 3 ? atoi(argv[3]) : 0; + int reps = argc > 4 ? atoi(argv[4]) : 3; + + bench_regions_t regs; + if (read_bed(bed, ®s) < 0) return 2; + fprintf(stderr, "loaded %d regions from %s\n", regs.n, bed); + + bench_sample_t scratch; + for (int w = 0; w < warmup; w++) { + if (run_single(fn, ®s, &scratch, 0) < 0) { fprintf(stderr, "warmup single failed\n"); return 2; } + if (run_multi (fn, ®s, &scratch, 0) < 0) { fprintf(stderr, "warmup multi failed\n"); return 2; } + } + + /* Untimed parity pass. Each path streams every emitted line through a + * running fnv1a64 hash and a byte counter; we then compare {length, + * hash} across paths. Constant memory regardless of panel size, so it + * scales to the 100k+ region case where holding two full byte streams + * (~1.7 GB each on a 1kGP VCF) would trigger swap or OOM. Doubles as + * a real warmup pass for the OS page cache, so the timed reps that + * follow operate on equivalent cache state. */ + printf("=== parity check ===\n"); + bench_sample_t parity_a, parity_b; + if (run_single(fn, ®s, &parity_a, 1) < 0) { + fprintf(stderr, "parity single failed\n"); return 2; + } + if (run_multi (fn, ®s, &parity_b, 1) < 0) { + fprintf(stderr, "parity multi failed\n"); return 2; + } + printf(" single: %ld records, %zu bytes, fnv1a=0x%016llx\n", + parity_a.records, parity_a.parity_bytes, + (unsigned long long) parity_a.parity_hash); + printf(" multi : %ld records, %zu bytes, fnv1a=0x%016llx\n", + parity_b.records, parity_b.parity_bytes, + (unsigned long long) parity_b.parity_hash); + + int rc = 0; + if (parity_a.parity_bytes != parity_b.parity_bytes || + parity_a.parity_hash != parity_b.parity_hash || + parity_a.records != parity_b.records) { + fprintf(stderr, "PARITY MISMATCH:" + " single=%ld records/%zu bytes/fnv1a=0x%016llx," + " multi=%ld records/%zu bytes/fnv1a=0x%016llx\n", + parity_a.records, parity_a.parity_bytes, + (unsigned long long) parity_a.parity_hash, + parity_b.records, parity_b.parity_bytes, + (unsigned long long) parity_b.parity_hash); + rc = 3; + } else { + printf(" PARITY OK: %ld records / %zu bytes match across paths\n", + parity_a.records, parity_a.parity_bytes); + } + if (rc) { free_regions(®s); return rc; } + + double sum_wall_single = 0, sum_wall_multi = 0; + long rec_single = parity_a.records, rec_multi = parity_b.records; + printf("=== bench_tbx_regions: %s with %d regions ===\n", fn, regs.n); + for (int r = 0; r < reps; r++) { + bench_sample_t a, b; + if (run_single(fn, ®s, &a, 0) < 0) return 2; + if (run_multi (fn, ®s, &b, 0) < 0) return 2; + printf("rep %d:\n", r); + print_sample("single", a); + print_sample("multi ", b); + printf(" speedup wall = %.2fx\n", a.wall_s / b.wall_s); + sum_wall_single += a.wall_s; + sum_wall_multi += b.wall_s; + if (a.records != rec_single || b.records != rec_multi) { + fprintf(stderr, "WARN: record count differs across reps " + "(single %ld -> %ld, multi %ld -> %ld)\n", + rec_single, a.records, rec_multi, b.records); + } + } + printf("=== summary ===\n"); + printf(" single avg wall = %.3fs over %d reps\n", sum_wall_single / reps, reps); + printf(" multi avg wall = %.3fs over %d reps\n", sum_wall_multi / reps, reps); + printf(" avg speedup = %.2fx\n", sum_wall_single / sum_wall_multi); + printf(" records single = %ld\n", rec_single); + printf(" records multi = %ld\n", rec_multi); + + free_regions(®s); + return 0; +} diff --git a/test/tbx-multi.bed b/test/tbx-multi.bed new file mode 100644 index 000000000..7e020468c --- /dev/null +++ b/test/tbx-multi.bed @@ -0,0 +1,6 @@ +1 99 100 region_a +1 149 150 region_b +1 199 200 region_c +1 299 300 region_d +2 99 100 region_e +2 199 200 region_f diff --git a/test/tbx-multi.gff b/test/tbx-multi.gff new file mode 100644 index 000000000..e1a20af6b --- /dev/null +++ b/test/tbx-multi.gff @@ -0,0 +1,7 @@ +##gff-version 3 +1 test feature 100 100 . + . ID=region_a +1 test feature 150 150 . + . ID=region_b +1 test feature 200 200 . + . ID=region_c +1 test feature 300 300 . + . ID=region_d +2 test feature 100 100 . + . ID=region_e +2 test feature 200 200 . + . ID=region_f diff --git a/test/tbx-multi.vcf b/test/tbx-multi.vcf new file mode 100644 index 000000000..44c4840b8 --- /dev/null +++ b/test/tbx-multi.vcf @@ -0,0 +1,10 @@ +##fileformat=VCFv4.2 +##contig= +##contig= +#CHROM POS ID REF ALT QUAL FILTER INFO +1 100 . A G 100 PASS . +1 150 . A G 100 PASS . +1 200 . A G 100 PASS . +1 300 . A G 100 PASS . +2 100 . A G 100 PASS . +2 200 . A G 100 PASS . diff --git a/test/test.pl b/test/test.pl index eaa65ea30..c5820d45f 100755 --- a/test/test.pl +++ b/test/test.pl @@ -59,6 +59,7 @@ run_test('test_bcf_sr_no_index',$opts); run_test('test_bcf_sr_range', $opts); run_test('test_bcf_sr_hreader', $opts); +run_test('test_tbx_multi', $opts); run_test('test_command',$opts,cmd=>'test-bcf-translate -',out=>'test-bcf-translate.out'); run_test('test_convert_padded_header',$opts); run_test('test_rebgzip',$opts); @@ -1425,6 +1426,35 @@ sub test_bcf_sr_hreader { passed($opts, $test); } +sub test_tbx_multi { + # Parity test for tbx_itr_regions (samtools/htslib PR for multi-region + # tabix iteration). Runs against three tabix preset formats (BED, VCF, + # GFF) to verify that the single-region path and the multi-region path + # produce byte-identical output across all of them. Plus a failure-path + # check for NULL inputs (ownership invariant). + my ($opts) = @_; + my $test = "test_tbx_multi"; + + my @fixtures = ( + { src => 'tbx-multi.bed', preset => 'bed' }, + { src => 'tbx-multi.vcf', preset => 'vcf' }, + { src => 'tbx-multi.gff', preset => 'gff' }, + ); + + foreach my $f (@fixtures) { + my $src = "$$opts{path}/$$f{src}"; + my $dst = "$$opts{tmp}/$$f{src}.gz"; + + my ($r1) = _cmd("$$opts{bin}/bgzip -c $src > $dst"); + if ($r1) { failed($opts, $test, "bgzip failed for $$f{src}: $r1"); return; } + my ($r2) = _cmd("$$opts{bin}/tabix -p $$f{preset} $dst"); + if ($r2) { failed($opts, $test, "tabix -p $$f{preset} failed for $$f{src}: $r2"); return; } + my ($r3) = _cmd("$$opts{path}/test_tbx_multi $dst"); + if ($r3) { failed($opts, $test, "test_tbx_multi exited $r3 for $$f{src} (parity mismatch?)"); return; } + } + passed($opts, $test); +} + sub test_command { my ($opts, %args) = @_; diff --git a/test/test_tbx_multi.c b/test/test_tbx_multi.c new file mode 100644 index 000000000..348dcc98e --- /dev/null +++ b/test/test_tbx_multi.c @@ -0,0 +1,393 @@ +/* test_tbx_multi.c -- correctness test for tbx_itr_regions. + + Compares the records returned by the multi-region path (tbx_itr_regions + plus hts_itr_multi_next) against the records returned by walking the + equivalent regions one at a time through the single-region path + (tbx_itr_queryi plus tbx_itr_next plus tbx_readrec). + + If both paths return byte-identical line sets, the new function is + correct by construction relative to the well-trodden single-region + callback. Otherwise the test reports the difference and exits non-zero. + + Copyright (C) 2026 Genome Research Ltd. + + Permission is hereby granted, free of charge, to any person obtaining a copy + of this software and associated documentation files (the "Software"), to deal + in the Software without restriction, including without limitation the rights + to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included in + all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + DEALINGS IN THE SOFTWARE. +*/ + +#include + +#include +#include +#include + +#include "../htslib/hts.h" +#include "../htslib/tbx.h" +#include "../htslib/kstring.h" + +/* Regions queried by both paths. Designed to exercise the multi-region code + * path under several shapes: point queries, wide range hitting multiple + * records, empty region returning nothing, and cross-contig spread. + * The fixture contains: + * 1\t99\t100\tregion_a (1-based pos 100) + * 1\t149\t150\tregion_b (pos 150) + * 1\t199\t200\tregion_c (pos 200) + * 1\t299\t300\tregion_d (pos 300) + * 2\t99\t100\tregion_e + * 2\t199\t200\tregion_f + */ +static const struct { const char *chrom; hts_pos_t beg, end; } QUERIES[] = { + { "1", 99, 100 }, /* point query, hits region_a */ + { "1", 140, 210 }, /* wide range, hits region_b and region_c */ + { "1", 5000, 6000 }, /* empty region, no records */ + { "1", 290, 310 }, /* point-ish query, hits region_d */ + { "2", 0, 1000 }, /* spans full contig 2 range, hits region_e and region_f */ +}; +static const int N_QUERIES = sizeof(QUERIES) / sizeof(QUERIES[0]); + +/* Append a line to a kstring. */ +static void append_line(kstring_t *out, const char *line) +{ + kputs(line, out); + kputc('\n', out); +} + +/* Single-region path: one tbx_itr_queryi per query, append results. */ +static int run_single_region(const char *filename, kstring_t *out) +{ + htsFile *fp = hts_open(filename, "r"); + if (!fp) { fprintf(stderr, "single: hts_open failed: %s\n", filename); return -1; } + tbx_t *tbx = tbx_index_load(filename); + if (!tbx) { fprintf(stderr, "single: tbx_index_load failed\n"); hts_close(fp); return -1; } + + kstring_t s = {0, 0, 0}; + int rc = 0; + for (int i = 0; i < N_QUERIES; i++) { + int tid = tbx_name2id(tbx, QUERIES[i].chrom); + if (tid < 0) { + fprintf(stderr, "single: no contig %s in fixture\n", QUERIES[i].chrom); + rc = -1; break; + } + hts_itr_t *itr = tbx_itr_queryi(tbx, tid, QUERIES[i].beg, QUERIES[i].end); + if (!itr) { + fprintf(stderr, "single: tbx_itr_queryi failed for %s:%lld-%lld\n", + QUERIES[i].chrom, (long long)QUERIES[i].beg, (long long)QUERIES[i].end); + rc = -1; break; + } + while (tbx_itr_next(fp, tbx, itr, &s) >= 0) { + append_line(out, s.s); + } + tbx_itr_destroy(itr); + } + free(s.s); + tbx_destroy(tbx); + hts_close(fp); + return rc; +} + +/* Multi-region path: build one reglist, one iterator via tbx_itr_regions, + * drain it. */ +static int run_multi_region(const char *filename, kstring_t *out) +{ + htsFile *fp = hts_open(filename, "r"); + if (!fp) { fprintf(stderr, "multi: hts_open failed: %s\n", filename); return -1; } + tbx_t *tbx = tbx_index_load(filename); + if (!tbx) { fprintf(stderr, "multi: tbx_index_load failed\n"); hts_close(fp); return -1; } + + /* Build a reglist: group queries by contig. hts_reglist_t holds intervals + * for ONE tid each, so one reglist per distinct chrom. */ + hts_reglist_t *reglist = calloc(N_QUERIES, sizeof(hts_reglist_t)); + if (!reglist) { + tbx_destroy(tbx); hts_close(fp); return -1; + } + int n_reg = 0; + for (int i = 0; i < N_QUERIES; i++) { + int tid = tbx_name2id(tbx, QUERIES[i].chrom); + if (tid < 0) { + fprintf(stderr, "multi: no contig %s in fixture\n", QUERIES[i].chrom); + free(reglist); tbx_destroy(tbx); hts_close(fp); return -1; + } + int idx = -1; + for (int j = 0; j < n_reg; j++) { + if (reglist[j].tid == tid) { idx = j; break; } + } + if (idx < 0) { + idx = n_reg++; + reglist[idx].tid = tid; + reglist[idx].reg = NULL; + reglist[idx].count = 0; + reglist[idx].intervals = NULL; + reglist[idx].min_beg = QUERIES[i].beg; + reglist[idx].max_end = QUERIES[i].end; + } + /* realloc into a temp so the original buffer is not lost on NULL. */ + hts_pair_pos_t *grown = realloc(reglist[idx].intervals, + (reglist[idx].count + 1) * sizeof(hts_pair_pos_t)); + if (!grown) { + for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + free(reglist); tbx_destroy(tbx); hts_close(fp); return -1; + } + reglist[idx].intervals = grown; + reglist[idx].intervals[reglist[idx].count].beg = QUERIES[i].beg; + reglist[idx].intervals[reglist[idx].count].end = QUERIES[i].end; + reglist[idx].count++; + if (QUERIES[i].beg < reglist[idx].min_beg) reglist[idx].min_beg = QUERIES[i].beg; + if (QUERIES[i].end > reglist[idx].max_end) reglist[idx].max_end = QUERIES[i].end; + } + + hts_itr_t *itr = tbx_itr_regions(fp, tbx, reglist, n_reg); + if (!itr) { + fprintf(stderr, "multi: tbx_itr_regions returned NULL\n"); + for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + free(reglist); + tbx_destroy(tbx); hts_close(fp); return -1; + } + + kstring_t s = {0, 0, 0}; + while (hts_itr_multi_next(fp, itr, &s) >= 0) { + append_line(out, s.s); + } + + free(s.s); + hts_itr_destroy(itr); + /* Caller always owns the original reglist (tbx_itr_regions deep copies). */ + for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + free(reglist); + tbx_destroy(tbx); + hts_close(fp); + return 0; +} + +/* Dedup test: tbx_itr_regions documents that overlapping and duplicate + * intervals yield each matching record at most once. Verify by constructing + * a reglist with three overlapping intervals that together cover the same + * three records four times over (each record is in two intervals at least), + * and asserting we still get exactly three records back. */ +static int test_dedup_and_overlap(const char *filename) +{ + htsFile *fp = hts_open(filename, "r"); + if (!fp) return -1; + tbx_t *tbx = tbx_index_load(filename); + if (!tbx) { hts_close(fp); return -1; } + + int tid = tbx_name2id(tbx, "1"); + if (tid < 0) { tbx_destroy(tbx); hts_close(fp); return -1; } + + /* Three intervals on contig 1, deliberately overlapping. The fixture + * has records at positions 100, 150, 200, 300 (1-based) on contig 1. + * [99, 200) covers records a (100), b (150), c (200-boundary) + * [120, 250) covers b, c + * [99, 100) covers a (full duplicate of the start of interval 1) + * If dedup is working, we expect exactly 3 records: a, b, c. */ + hts_reglist_t *reg = calloc(1, sizeof(hts_reglist_t)); + if (!reg) { tbx_destroy(tbx); hts_close(fp); return -1; } + reg[0].tid = tid; + reg[0].count = 3; + reg[0].intervals = calloc(3, sizeof(hts_pair_pos_t)); + if (!reg[0].intervals) { free(reg); tbx_destroy(tbx); hts_close(fp); return -1; } + reg[0].intervals[0].beg = 99; reg[0].intervals[0].end = 200; + reg[0].intervals[1].beg = 120; reg[0].intervals[1].end = 250; + reg[0].intervals[2].beg = 99; reg[0].intervals[2].end = 100; + reg[0].min_beg = 99; reg[0].max_end = 250; + + hts_itr_t *itr = tbx_itr_regions(fp, tbx, reg, 1); + if (!itr) { + fprintf(stderr, "FAIL: dedup test could not construct iterator\n"); + free(reg[0].intervals); free(reg); + tbx_destroy(tbx); hts_close(fp); + return -1; + } + + kstring_t s = {0, 0, 0}; + int n = 0; + while (hts_itr_multi_next(fp, itr, &s) >= 0) n++; + free(s.s); + hts_itr_destroy(itr); + /* Caller always owns the original reglist (tbx_itr_regions deep copies). */ + free(reg[0].intervals); free(reg); + tbx_destroy(tbx); + hts_close(fp); + + /* Three records expected: a, b, c. If dedup were broken we'd get 5+. */ + if (n != 3) { + fprintf(stderr, "FAIL: dedup expected 3 records, got %d\n", n); + return -1; + } + fprintf(stderr, "OK: dedup yields exactly 3 records across overlapping intervals\n"); + return 0; +} + +/* Failure-path test: confirm tbx_itr_regions returns NULL cleanly when given + * bad inputs. + * + * Per tbx.h, the caller ALWAYS owns the original reglist regardless of + * outcome (the function deep copies internally). This test exercises the + * input-validation guard plus a successful-then-freed case, and the + * ASan-instrumented build will catch any double-free or use-after-free. */ +static int test_failure_paths(const char *filename) +{ + htsFile *fp = hts_open(filename, "r"); + if (!fp) return -1; + tbx_t *tbx = tbx_index_load(filename); + if (!tbx) { hts_close(fp); return -1; } + + hts_reglist_t *reg = calloc(1, sizeof(hts_reglist_t)); + if (!reg) { tbx_destroy(tbx); hts_close(fp); return -1; } + reg[0].tid = tbx_name2id(tbx, "1"); + reg[0].count = 1; + reg[0].intervals = calloc(1, sizeof(hts_pair_pos_t)); + if (!reg[0].intervals) { free(reg); tbx_destroy(tbx); hts_close(fp); return -1; } + reg[0].intervals[0].beg = 99; + reg[0].intervals[0].end = 100; + reg[0].min_beg = 99; + reg[0].max_end = 100; + + int rc = 0; + hts_itr_t *itr; + + struct { const char *name; htsFile *fp; tbx_t *tbx; hts_reglist_t *reg; int count; } cases[] = { + { "NULL htsFile", NULL, tbx, reg, 1 }, + { "NULL tbx", fp, NULL, reg, 1 }, + { "NULL fp+tbx", NULL, NULL, reg, 1 }, + { "NULL reglist", fp, tbx, NULL, 1 }, + { "count == 0", fp, tbx, reg, 0 }, + { "negative count", fp, tbx, reg, -1 }, + }; + const int n_cases = sizeof(cases) / sizeof(cases[0]); + + for (int i = 0; i < n_cases; i++) { + itr = tbx_itr_regions(cases[i].fp, cases[i].tbx, cases[i].reg, cases[i].count); + if (itr != NULL) { + fprintf(stderr, "FAIL: %s should return NULL, got non-NULL\n", cases[i].name); + hts_itr_destroy(itr); + rc = 1; + } else { + fprintf(stderr, "OK: %s returns NULL\n", cases[i].name); + } + } + + /* Test the .reg name-string resolution paths added by the deep-copy + * refactor. Both should succeed in constructing an iterator: one + * resolves to a known tid, one warns and resolves to a negative tid + * (the iterator still constructs but yields no records for that + * entry). Caller still owns reg afterwards. */ + { + hts_reglist_t *named = calloc(1, sizeof(*named)); + named[0].reg = "1"; + named[0].count = 1; + named[0].intervals = calloc(1, sizeof(hts_pair_pos_t)); + named[0].intervals[0].beg = 99; + named[0].intervals[0].end = 100; + named[0].min_beg = 99; named[0].max_end = 100; + + itr = tbx_itr_regions(fp, tbx, named, 1); + if (!itr) { + fprintf(stderr, "FAIL: .reg=\"1\" should resolve and succeed\n"); + rc = 1; + } else { + fprintf(stderr, "OK: .reg=\"1\" resolved via tbx_name2id\n"); + hts_itr_destroy(itr); + } + free(named[0].intervals); free(named); + } + { + hts_reglist_t *named = calloc(1, sizeof(*named)); + named[0].reg = "nonexistent_contig"; + named[0].count = 1; + named[0].intervals = calloc(1, sizeof(hts_pair_pos_t)); + named[0].intervals[0].beg = 99; + named[0].intervals[0].end = 100; + named[0].min_beg = 99; named[0].max_end = 100; + + itr = tbx_itr_regions(fp, tbx, named, 1); + if (!itr) { + fprintf(stderr, "OK: unknown .reg name returns NULL (acceptable)\n"); + } else { + fprintf(stderr, "OK: unknown .reg name resolved to negative tid" + " and iterator constructed (will yield nothing)\n"); + hts_itr_destroy(itr); + } + free(named[0].intervals); free(named); + } + + /* The ASan-instrumented build catches any double-free or + * use-after-free; the deep-copy contract means this trailing free is + * always safe regardless of which paths above hit. */ + free(reg[0].intervals); + free(reg); + + tbx_destroy(tbx); + hts_close(fp); + return rc; +} + +int main(int argc, char **argv) +{ + if (argc != 2) { + fprintf(stderr, "Usage: %s file.tabix.gz\n", argv[0]); + return 1; + } + const char *filename = argv[1]; + + kstring_t single_out = {0, 0, 0}; + kstring_t multi_out = {0, 0, 0}; + + if (run_single_region(filename, &single_out) < 0) { + fprintf(stderr, "single-region run failed\n"); + free(single_out.s); free(multi_out.s); + return 2; + } + if (run_multi_region(filename, &multi_out) < 0) { + fprintf(stderr, "multi-region run failed\n"); + free(single_out.s); free(multi_out.s); + return 2; + } + + fprintf(stderr, "single-region output (%zu bytes):\n%s", single_out.l, single_out.s); + fprintf(stderr, "multi-region output (%zu bytes):\n%s", multi_out.l, multi_out.s); + + int rc = 0; + if (single_out.l != multi_out.l || memcmp(single_out.s, multi_out.s, single_out.l) != 0) { + fprintf(stderr, "PARITY MISMATCH between single-region and multi-region outputs\n"); + rc = 1; + } else { + fprintf(stderr, "PARITY OK: %zu bytes identical\n", single_out.l); + } + + /* Dedup coverage (only for BED fixture, which has the positions the + * test expects; VCF/GFF fixtures share positions so it works for them + * too, but the constants in the test are BED-centric). */ + if (test_dedup_and_overlap(filename) != 0) { + fprintf(stderr, "DEDUP/OVERLAP TESTS failed\n"); + rc = 1; + } else { + fprintf(stderr, "DEDUP/OVERLAP TESTS OK\n"); + } + + /* Failure-path coverage */ + if (test_failure_paths(filename) != 0) { + fprintf(stderr, "FAILURE-PATH TESTS failed\n"); + rc = 1; + } else { + fprintf(stderr, "FAILURE-PATH TESTS OK\n"); + } + + free(single_out.s); + free(multi_out.s); + return rc; +} From 71d36834e89d38891807d9ec7153ebd87acbecef Mon Sep 17 00:00:00 2001 From: Carsten Erickson Date: Wed, 27 May 2026 08:27:05 -0700 Subject: [PATCH 2/2] Fix CI: gnu90 compatibility and gitignore for built test binaries Cirrus reported four configurations failing on the PR head: rockylinux-gcc builds with -std=gnu90 -Werror and rejected the for-loop initial declarations and mid-block variable declarations introduced in tbx.c, test_tbx_multi.c, and bench_tbx_regions.c. debian-gcc, ubuntu-arm, ubuntu-clang three configs run a DO_UNTRACKED_FILE_CHECK that flagged the two new test binaries (test/test_tbx_multi and test/bench_tbx_regions) as untracked. This commit: - Moves all variable declarations to block heads in tbx.c (in tbx_reglist_dup and tbx_itr_regions), in test/test_tbx_multi.c, and in test/bench_tbx_regions.c, so the files compile under -std=gnu90. - Adds /test/bench_tbx_regions and /test/test_tbx_multi entries to .gitignore, alphabetised next to the existing test binary entries. Verified locally: gcc -std=gnu90 -Wall -Wformat=2 -Wextra -Werror on all three .c files reports no errors, and the full ASan + UBSan htslib test suite (361 tests) passes. Signed-off-by: Carsten Erickson Assisted-by: Claude (claude-opus-4-7) --- .gitignore | 2 + tbx.c | 34 +++++++----- test/bench_tbx_regions.c | 108 ++++++++++++++++++++++--------------- test/test_tbx_multi.c | 111 +++++++++++++++++++++++---------------- 4 files changed, 155 insertions(+), 100 deletions(-) diff --git a/.gitignore b/.gitignore index a485dce54..e88c16aaf 100644 --- a/.gitignore +++ b/.gitignore @@ -59,6 +59,7 @@ shlib-exports-*.txt /test/sam /test/ref_cache/*.tmp.* /test/tabix/*.tmp.* +/test/bench_tbx_regions /test/test_alloc /test/test-bcf-sr /test/test-bcf-translate @@ -78,6 +79,7 @@ shlib-exports-*.txt /test/test_realn /test/test-regidx /test/test_str2int +/test/test_tbx_multi /test/test_time_funcs /test/test-vcf-api /test/test-vcf-sweep diff --git a/tbx.c b/tbx.c index 61a50a2ab..bbe29beac 100644 --- a/tbx.c +++ b/tbx.c @@ -439,10 +439,16 @@ static int tbx_readrec_multi(BGZF *fp, void *ignored, void *sv, static int tbx_reglist_dup(tbx_t *tbx, const hts_reglist_t *src, int count, hts_reglist_t **out) { - hts_reglist_t *dst = calloc(count, sizeof(*dst)); + hts_reglist_t *dst; + int i; + + dst = calloc(count, sizeof(*dst)); if (!dst) return -1; - for (int i = 0; i < count; i++) { + for (i = 0; i < count; i++) { + uint32_t n, k; + hts_pos_t mn, mx; + dst[i].tid = src[i].tid; dst[i].count = src[i].count; dst[i].min_beg = src[i].min_beg; @@ -476,7 +482,7 @@ static int tbx_reglist_dup(tbx_t *tbx, const hts_reglist_t *src, int count, * mismatch between the size we allocate and the size we believe * dst holds, even though callers are not expected to mutate src * concurrently. */ - uint32_t n = src[i].count; + n = src[i].count; dst[i].intervals = malloc((size_t) n * sizeof(hts_pair_pos_t)); if (!dst[i].intervals) { hts_reglist_free(dst, count); @@ -489,9 +495,9 @@ static int tbx_reglist_dup(tbx_t *tbx, const hts_reglist_t *src, int count, /* Recompute min_beg / max_end from the copied intervals so they * remain consistent even if the caller's summary fields were * stale or if the later qsort reorders intervals[0]. */ - hts_pos_t mn = dst[i].intervals[0].beg; - hts_pos_t mx = dst[i].intervals[0].end; - for (uint32_t k = 1; k < n; k++) { + mn = dst[i].intervals[0].beg; + mx = dst[i].intervals[0].end; + for (k = 1; k < n; k++) { if (dst[i].intervals[k].beg < mn) mn = dst[i].intervals[k].beg; if (dst[i].intervals[k].end > mx) mx = dst[i].intervals[k].end; } @@ -506,15 +512,19 @@ static int tbx_reglist_dup(tbx_t *tbx, const hts_reglist_t *src, int count, hts_itr_t *tbx_itr_regions(htsFile *fp, tbx_t *tbx, hts_reglist_t *reglist, int count) { + BGZF *bgzf; + hts_reglist_t *owned = NULL; + hts_itr_t *itr; + int i; + if (!fp || !tbx || !reglist || count <= 0) return NULL; - BGZF *bgzf = hts_get_bgzfp(fp); + bgzf = hts_get_bgzfp(fp); if (!bgzf) return NULL; /* Make our own copy of the reglist so the caller's input is never * mutated and never assumed to be live past this call. The copy also * lets us pre-resolve any .reg name strings to tids so we can pass * getid=NULL safely to hts_itr_regions below. */ - hts_reglist_t *owned = NULL; if (tbx_reglist_dup(tbx, reglist, count, &owned) < 0) return NULL; @@ -522,16 +532,16 @@ hts_itr_t *tbx_itr_regions(htsFile *fp, tbx_t *tbx, * within each reglist entry to be sorted by start position; its * advance logic uses the interval index packed into the offset list * which relies on sorted iteration. Sort our copy now. */ - for (int i = 0; i < count; i++) { + for (i = 0; i < count; i++) { if (owned[i].intervals && owned[i].count > 1) { qsort(owned[i].intervals, owned[i].count, sizeof(hts_pair_pos_t), compare_hts_pair_pos_t); } } - hts_itr_t *itr = hts_itr_regions(tbx->idx, owned, count, NULL, NULL, - hts_itr_multi_bam, tbx_readrec_multi, - tbx_pseek, tbx_ptell); + itr = hts_itr_regions(tbx->idx, owned, count, NULL, NULL, + hts_itr_multi_bam, tbx_readrec_multi, + tbx_pseek, tbx_ptell); if (!itr) { /* hts_itr_regions has exactly one failure path that does not * free reglist via hts_itr_destroy: the OOM at its initial diff --git a/test/bench_tbx_regions.c b/test/bench_tbx_regions.c index d60dc5bc0..a617f6260 100644 --- a/test/bench_tbx_regions.c +++ b/test/bench_tbx_regions.c @@ -65,7 +65,8 @@ static inline uint64_t fnv1a64_update(uint64_t h, const void *data, size_t len) { const unsigned char *p = (const unsigned char *) data; - for (size_t i = 0; i < len; i++) { + size_t i; + for (i = 0; i < len; i++) { h ^= p[i]; h *= FNV1A64_PRIME; } @@ -134,7 +135,8 @@ static int read_bed(const char *path, bench_regions_t *out) static void free_regions(bench_regions_t *r) { - for (int i = 0; i < r->n; i++) free(r->items[i].chrom); + int i; + for (i = 0; i < r->n; i++) free(r->items[i].chrom); free(r->items); } @@ -147,10 +149,8 @@ static void free_regions(bench_regions_t *r) static int run_single(const char *fn, const bench_regions_t *regs, bench_sample_t *s, int with_parity) { - htsFile *fp = hts_open(fn, "r"); - if (!fp) return -1; - tbx_t *tbx = tbx_index_load(fn); - if (!tbx) { hts_close(fp); return -1; } + htsFile *fp; + tbx_t *tbx; kstring_t line = {0,0,0}; long records = 0; uint64_t hash = FNV1A64_OFFSET; @@ -158,13 +158,21 @@ static int run_single(const char *fn, const bench_regions_t *regs, static const char newline = '\n'; struct rusage ru0, ru1; struct timespec t0, t1; + int i; + + fp = hts_open(fn, "r"); + if (!fp) return -1; + tbx = tbx_index_load(fn); + if (!tbx) { hts_close(fp); return -1; } getrusage(RUSAGE_SELF, &ru0); clock_gettime(CLOCK_MONOTONIC, &t0); - for (int i = 0; i < regs->n; i++) { - int tid = tbx_name2id(tbx, regs->items[i].chrom); + for (i = 0; i < regs->n; i++) { + int tid; + hts_itr_t *itr; + tid = tbx_name2id(tbx, regs->items[i].chrom); if (tid < 0) continue; - hts_itr_t *itr = tbx_itr_queryi(tbx, tid, regs->items[i].beg, - regs->items[i].end); + itr = tbx_itr_queryi(tbx, tid, regs->items[i].beg, + regs->items[i].end); if (!itr) continue; while (tbx_itr_next(fp, tbx, itr, &line) >= 0) { records++; @@ -195,20 +203,35 @@ static int run_single(const char *fn, const bench_regions_t *regs, static int run_multi(const char *fn, const bench_regions_t *regs, bench_sample_t *s, int with_parity) { - htsFile *fp = hts_open(fn, "r"); + htsFile *fp; + tbx_t *tbx; + hts_reglist_t *reglist; + hts_itr_t *itr; + kstring_t line = {0,0,0}; + long records = 0; + uint64_t hash = FNV1A64_OFFSET; + size_t bytes = 0; + static const char newline = '\n'; + struct rusage ru0, ru1; + struct timespec t0, t1; + int n_reg = 0; + int i, j; + + fp = hts_open(fn, "r"); if (!fp) return -1; - tbx_t *tbx = tbx_index_load(fn); + tbx = tbx_index_load(fn); if (!tbx) { hts_close(fp); return -1; } /* Group regions by tid into one reglist entry per contig. */ - hts_reglist_t *reglist = calloc(regs->n, sizeof(*reglist)); + reglist = calloc(regs->n, sizeof(*reglist)); if (!reglist) { tbx_destroy(tbx); hts_close(fp); return -1; } - int n_reg = 0; - for (int i = 0; i < regs->n; i++) { - int tid = tbx_name2id(tbx, regs->items[i].chrom); + for (i = 0; i < regs->n; i++) { + int tid, idx; + hts_pair_pos_t *grown; + tid = tbx_name2id(tbx, regs->items[i].chrom); if (tid < 0) continue; - int idx = -1; - for (int j = 0; j < n_reg; j++) { + idx = -1; + for (j = 0; j < n_reg; j++) { if (reglist[j].tid == tid) { idx = j; break; } } if (idx < 0) { @@ -220,10 +243,10 @@ static int run_multi(const char *fn, const bench_regions_t *regs, reglist[idx].min_beg = regs->items[i].beg; reglist[idx].max_end = regs->items[i].end; } - hts_pair_pos_t *grown = realloc(reglist[idx].intervals, - (reglist[idx].count + 1) * sizeof(hts_pair_pos_t)); + grown = realloc(reglist[idx].intervals, + (reglist[idx].count + 1) * sizeof(hts_pair_pos_t)); if (!grown) { - for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + for (j = 0; j < n_reg; j++) free(reglist[j].intervals); free(reglist); tbx_destroy(tbx); hts_close(fp); return -1; } reglist[idx].intervals = grown; @@ -236,19 +259,12 @@ static int run_multi(const char *fn, const bench_regions_t *regs, reglist[idx].max_end = regs->items[i].end; } - kstring_t line = {0,0,0}; - long records = 0; - uint64_t hash = FNV1A64_OFFSET; - size_t bytes = 0; - static const char newline = '\n'; - struct rusage ru0, ru1; - struct timespec t0, t1; getrusage(RUSAGE_SELF, &ru0); clock_gettime(CLOCK_MONOTONIC, &t0); - hts_itr_t *itr = tbx_itr_regions(fp, tbx, reglist, n_reg); + itr = tbx_itr_regions(fp, tbx, reglist, n_reg); if (!itr) { clock_gettime(CLOCK_MONOTONIC, &t1); - for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + for (j = 0; j < n_reg; j++) free(reglist[j].intervals); free(reglist); free(line.s); tbx_destroy(tbx); hts_close(fp); return -1; @@ -284,21 +300,31 @@ static void print_sample(const char *label, bench_sample_t s) int main(int argc, char **argv) { + const char *fn; + const char *bed; + int warmup; + int reps; + bench_regions_t regs; + bench_sample_t scratch; + bench_sample_t parity_a, parity_b; + double sum_wall_single = 0, sum_wall_multi = 0; + long rec_single, rec_multi; + int rc = 0; + int w, r; + if (argc < 3) { fprintf(stderr, "Usage: %s FILE.gz REGIONS.bed [WARMUP] [REPS]\n", argv[0]); return 1; } - const char *fn = argv[1]; - const char *bed = argv[2]; - int warmup = argc > 3 ? atoi(argv[3]) : 0; - int reps = argc > 4 ? atoi(argv[4]) : 3; + fn = argv[1]; + bed = argv[2]; + warmup = argc > 3 ? atoi(argv[3]) : 0; + reps = argc > 4 ? atoi(argv[4]) : 3; - bench_regions_t regs; if (read_bed(bed, ®s) < 0) return 2; fprintf(stderr, "loaded %d regions from %s\n", regs.n, bed); - bench_sample_t scratch; - for (int w = 0; w < warmup; w++) { + for (w = 0; w < warmup; w++) { if (run_single(fn, ®s, &scratch, 0) < 0) { fprintf(stderr, "warmup single failed\n"); return 2; } if (run_multi (fn, ®s, &scratch, 0) < 0) { fprintf(stderr, "warmup multi failed\n"); return 2; } } @@ -311,7 +337,6 @@ int main(int argc, char **argv) * a real warmup pass for the OS page cache, so the timed reps that * follow operate on equivalent cache state. */ printf("=== parity check ===\n"); - bench_sample_t parity_a, parity_b; if (run_single(fn, ®s, &parity_a, 1) < 0) { fprintf(stderr, "parity single failed\n"); return 2; } @@ -325,7 +350,6 @@ int main(int argc, char **argv) parity_b.records, parity_b.parity_bytes, (unsigned long long) parity_b.parity_hash); - int rc = 0; if (parity_a.parity_bytes != parity_b.parity_bytes || parity_a.parity_hash != parity_b.parity_hash || parity_a.records != parity_b.records) { @@ -343,10 +367,10 @@ int main(int argc, char **argv) } if (rc) { free_regions(®s); return rc; } - double sum_wall_single = 0, sum_wall_multi = 0; - long rec_single = parity_a.records, rec_multi = parity_b.records; + rec_single = parity_a.records; + rec_multi = parity_b.records; printf("=== bench_tbx_regions: %s with %d regions ===\n", fn, regs.n); - for (int r = 0; r < reps; r++) { + for (r = 0; r < reps; r++) { bench_sample_t a, b; if (run_single(fn, ®s, &a, 0) < 0) return 2; if (run_multi (fn, ®s, &b, 0) < 0) return 2; diff --git a/test/test_tbx_multi.c b/test/test_tbx_multi.c index 348dcc98e..be4835ba8 100644 --- a/test/test_tbx_multi.c +++ b/test/test_tbx_multi.c @@ -70,20 +70,26 @@ static void append_line(kstring_t *out, const char *line) /* Single-region path: one tbx_itr_queryi per query, append results. */ static int run_single_region(const char *filename, kstring_t *out) { - htsFile *fp = hts_open(filename, "r"); + htsFile *fp; + tbx_t *tbx; + kstring_t s = {0, 0, 0}; + int rc = 0; + int i; + + fp = hts_open(filename, "r"); if (!fp) { fprintf(stderr, "single: hts_open failed: %s\n", filename); return -1; } - tbx_t *tbx = tbx_index_load(filename); + tbx = tbx_index_load(filename); if (!tbx) { fprintf(stderr, "single: tbx_index_load failed\n"); hts_close(fp); return -1; } - kstring_t s = {0, 0, 0}; - int rc = 0; - for (int i = 0; i < N_QUERIES; i++) { - int tid = tbx_name2id(tbx, QUERIES[i].chrom); + for (i = 0; i < N_QUERIES; i++) { + int tid; + hts_itr_t *itr; + tid = tbx_name2id(tbx, QUERIES[i].chrom); if (tid < 0) { fprintf(stderr, "single: no contig %s in fixture\n", QUERIES[i].chrom); rc = -1; break; } - hts_itr_t *itr = tbx_itr_queryi(tbx, tid, QUERIES[i].beg, QUERIES[i].end); + itr = tbx_itr_queryi(tbx, tid, QUERIES[i].beg, QUERIES[i].end); if (!itr) { fprintf(stderr, "single: tbx_itr_queryi failed for %s:%lld-%lld\n", QUERIES[i].chrom, (long long)QUERIES[i].beg, (long long)QUERIES[i].end); @@ -104,26 +110,35 @@ static int run_single_region(const char *filename, kstring_t *out) * drain it. */ static int run_multi_region(const char *filename, kstring_t *out) { - htsFile *fp = hts_open(filename, "r"); + htsFile *fp; + tbx_t *tbx; + hts_reglist_t *reglist; + hts_itr_t *itr; + kstring_t s = {0, 0, 0}; + int n_reg = 0; + int i, j; + + fp = hts_open(filename, "r"); if (!fp) { fprintf(stderr, "multi: hts_open failed: %s\n", filename); return -1; } - tbx_t *tbx = tbx_index_load(filename); + tbx = tbx_index_load(filename); if (!tbx) { fprintf(stderr, "multi: tbx_index_load failed\n"); hts_close(fp); return -1; } /* Build a reglist: group queries by contig. hts_reglist_t holds intervals * for ONE tid each, so one reglist per distinct chrom. */ - hts_reglist_t *reglist = calloc(N_QUERIES, sizeof(hts_reglist_t)); + reglist = calloc(N_QUERIES, sizeof(hts_reglist_t)); if (!reglist) { tbx_destroy(tbx); hts_close(fp); return -1; } - int n_reg = 0; - for (int i = 0; i < N_QUERIES; i++) { - int tid = tbx_name2id(tbx, QUERIES[i].chrom); + for (i = 0; i < N_QUERIES; i++) { + int tid, idx; + hts_pair_pos_t *grown; + tid = tbx_name2id(tbx, QUERIES[i].chrom); if (tid < 0) { fprintf(stderr, "multi: no contig %s in fixture\n", QUERIES[i].chrom); free(reglist); tbx_destroy(tbx); hts_close(fp); return -1; } - int idx = -1; - for (int j = 0; j < n_reg; j++) { + idx = -1; + for (j = 0; j < n_reg; j++) { if (reglist[j].tid == tid) { idx = j; break; } } if (idx < 0) { @@ -136,10 +151,10 @@ static int run_multi_region(const char *filename, kstring_t *out) reglist[idx].max_end = QUERIES[i].end; } /* realloc into a temp so the original buffer is not lost on NULL. */ - hts_pair_pos_t *grown = realloc(reglist[idx].intervals, - (reglist[idx].count + 1) * sizeof(hts_pair_pos_t)); + grown = realloc(reglist[idx].intervals, + (reglist[idx].count + 1) * sizeof(hts_pair_pos_t)); if (!grown) { - for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + for (j = 0; j < n_reg; j++) free(reglist[j].intervals); free(reglist); tbx_destroy(tbx); hts_close(fp); return -1; } reglist[idx].intervals = grown; @@ -150,15 +165,14 @@ static int run_multi_region(const char *filename, kstring_t *out) if (QUERIES[i].end > reglist[idx].max_end) reglist[idx].max_end = QUERIES[i].end; } - hts_itr_t *itr = tbx_itr_regions(fp, tbx, reglist, n_reg); + itr = tbx_itr_regions(fp, tbx, reglist, n_reg); if (!itr) { fprintf(stderr, "multi: tbx_itr_regions returned NULL\n"); - for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + for (j = 0; j < n_reg; j++) free(reglist[j].intervals); free(reglist); tbx_destroy(tbx); hts_close(fp); return -1; } - kstring_t s = {0, 0, 0}; while (hts_itr_multi_next(fp, itr, &s) >= 0) { append_line(out, s.s); } @@ -166,7 +180,7 @@ static int run_multi_region(const char *filename, kstring_t *out) free(s.s); hts_itr_destroy(itr); /* Caller always owns the original reglist (tbx_itr_regions deep copies). */ - for (int j = 0; j < n_reg; j++) free(reglist[j].intervals); + for (j = 0; j < n_reg; j++) free(reglist[j].intervals); free(reglist); tbx_destroy(tbx); hts_close(fp); @@ -241,12 +255,19 @@ static int test_dedup_and_overlap(const char *filename) * ASan-instrumented build will catch any double-free or use-after-free. */ static int test_failure_paths(const char *filename) { - htsFile *fp = hts_open(filename, "r"); + htsFile *fp; + tbx_t *tbx; + hts_reglist_t *reg; + hts_itr_t *itr; + int rc = 0; + int i, n_cases; + + fp = hts_open(filename, "r"); if (!fp) return -1; - tbx_t *tbx = tbx_index_load(filename); + tbx = tbx_index_load(filename); if (!tbx) { hts_close(fp); return -1; } - hts_reglist_t *reg = calloc(1, sizeof(hts_reglist_t)); + reg = calloc(1, sizeof(hts_reglist_t)); if (!reg) { tbx_destroy(tbx); hts_close(fp); return -1; } reg[0].tid = tbx_name2id(tbx, "1"); reg[0].count = 1; @@ -257,27 +278,25 @@ static int test_failure_paths(const char *filename) reg[0].min_beg = 99; reg[0].max_end = 100; - int rc = 0; - hts_itr_t *itr; - - struct { const char *name; htsFile *fp; tbx_t *tbx; hts_reglist_t *reg; int count; } cases[] = { - { "NULL htsFile", NULL, tbx, reg, 1 }, - { "NULL tbx", fp, NULL, reg, 1 }, - { "NULL fp+tbx", NULL, NULL, reg, 1 }, - { "NULL reglist", fp, tbx, NULL, 1 }, - { "count == 0", fp, tbx, reg, 0 }, - { "negative count", fp, tbx, reg, -1 }, - }; - const int n_cases = sizeof(cases) / sizeof(cases[0]); - - for (int i = 0; i < n_cases; i++) { - itr = tbx_itr_regions(cases[i].fp, cases[i].tbx, cases[i].reg, cases[i].count); - if (itr != NULL) { - fprintf(stderr, "FAIL: %s should return NULL, got non-NULL\n", cases[i].name); - hts_itr_destroy(itr); - rc = 1; - } else { - fprintf(stderr, "OK: %s returns NULL\n", cases[i].name); + { + struct { const char *name; htsFile *fp; tbx_t *tbx; hts_reglist_t *reg; int count; } cases[] = { + { "NULL htsFile", NULL, tbx, reg, 1 }, + { "NULL tbx", fp, NULL, reg, 1 }, + { "NULL fp+tbx", NULL, NULL, reg, 1 }, + { "NULL reglist", fp, tbx, NULL, 1 }, + { "count == 0", fp, tbx, reg, 0 }, + { "negative count", fp, tbx, reg, -1 }, + }; + n_cases = sizeof(cases) / sizeof(cases[0]); + for (i = 0; i < n_cases; i++) { + itr = tbx_itr_regions(cases[i].fp, cases[i].tbx, cases[i].reg, cases[i].count); + if (itr != NULL) { + fprintf(stderr, "FAIL: %s should return NULL, got non-NULL\n", cases[i].name); + hts_itr_destroy(itr); + rc = 1; + } else { + fprintf(stderr, "OK: %s returns NULL\n", cases[i].name); + } } }