Add an inline version of bgzf_read for small reads. by jkbonfield · Pull Request #1772 · samtools/htslib

jkbonfield · 2024-04-19T11:37:16Z

The bgzf_read function is long and not something that's likely to get inlined even if it was in the header.

However most of the time our calls request a small amount of data and they fit within the buffer we've read, so we offer a static inline to do the memcpy when we can, falling back to the long function when we cannot.

In terms of CPU time it's not much difference, but the key thing is that it's often CPU time saved in a main thread given the bulk of the decode is often threaded. An example of test_view -B -@16

develop:

real    0m48.158s
user    6m2.901s
sys     0m28.134s

real    0m48.730s
user    6m3.707s
sys     0m28.473s

real    0m48.653s
user    6m5.215s
sys     0m28.637s

This PR:

real    0m41.731s
user    5m59.780s
sys     0m30.393s

real    0m41.945s
user    6m0.367s
sys     0m30.426s

So we can see it's a consistent win when threading, potentially 10-15% faster throughput depending on work loads.

jkbonfield · 2024-04-19T14:46:50Z

/usr/bin/time was reporting ~810% CPU utilisation for develop and ~940% for this PR, demonstrating the reduced time in main thread.

Commenting out the sanity checks in bam_read1 for if (c->n_cigar > 0) onwards increases that CPU utilisation to 1060% (1120% if I get rid of the tag2cigar call), leading to 37s real time. This shows that there is room for improvement, as these sort of checks ought to be executed in the multi-threaded decode instead. We thread the bgzf_read and not the bam_read.

This is why multi-threaded decoding of SAM and SAM.gz is sometimes more performant than BAM, despite the string->integer conversion overheads, as SAM parsing is entirely threadable.

daviesrob · 2024-05-09T09:00:42Z

-        for (i = 0; i < 8; ++i) ed_swap_4p(x + i);
+    if (fp->block_length - fp->block_offset > 32) {
+        // Avoid bgzf_read and a temporary copy to a local buffer
+        uint8_t *x = (uint8_t *)fp->uncompressed_block + fp->block_offset;


This shadows the earlier x, and also adds an almost-duplicate block to fill out the BAM structure. A way to avoid both problems might be something like this:

uint8_t tmp[32], *x; uint32_t new_l_data; // ... if (fp->block_length - fp->block_offset > 32) { x = (uint8_t *)fp->uncompressed_block + fp->block_offset; fp->block_offset += 32; } else { x = tmp; if (bgzf_read(fp, x, 32) != 32) return -3; } c->tid = le_to_u32(x); c->pos = le_to_i32(x+4); // ... etc.

I see what you mean by duplicate now. It's totally different code and not duplicated at all in style nor in implementation, but the end result is still the same thing. I've taken your suggestion and refactored it.

daviesrob · 2024-05-09T09:08:54Z

+        uint32_t x2   = le_to_u32(x+8);
+        c->bin        = x2>>16;
+        c->qual       = x2>>8&0xff;
+        c->l_qname    = x2&0xff;


The intermediate variable could be avoided by:

c->l_qname = le_to_u8(x + 8); c->qual = le_to_u8(x + 9); c->bin = le_to_u16(x + 10);

I'm not sure if it would be any quicker though...

I was mainly mirroring what it did before, but yes it may make things better (or worse). I'll experiment.

daviesrob · 2024-05-09T09:14:05Z

+        uint32_t x3   = le_to_u32(x+12);
+        c->flag       = x3>>16;
+        c->n_cigar    = x3&0xffff;


Similarly, could be:

c->n_cigar = le_to_u16(x + 12); c->flag = le_to_u16(x + 14);

daviesrob · 2024-05-09T09:34:27Z

 * Note a second interface that returns a bam pointer instead would avoid bam_copy1
 * in multi-threaded handling.  This may be worth considering for htslib2.
 */
+#define bgzf_read bgzf_read_small


I think it would be less obscure to replace the individual bgzf_read() calls, as there aren't many of them. Also bgzf_read(fp, x, 32) should remain as-is as the amount of data available would be known to be too small in that case.

You could also merge bgzf_read(fp, &block_len, 4) into the one for the rest of the core header. That would eliminate an entire bgzf_read call, and I don't really see a good reason why it wouldn't work as long as we were careful about checking error cases and returning the right values.

I'd previously tested merging the 4 byte and 32 byte read together and concluded it made no difference in speed (once I'd put the inline version in - it did do before that), and it meant more changes which would have made reviewing more work.

I can try it again if you want, but it's more changes to evaluate.

Also, I don't get the comment about bgzf_read(fp, x, 32) remaining as-is. Initially it wasn't because it's small enough to be inlined and a fixed 32-byte read may still cause the memcpy to be optimised away (I don't know), and because we weren't subverting things and copying out the struct directly initially. However since then, it may not matter much either way as that code path will be very rare. I think it's probably an irrelevance at best.

You also have to think about what's actually happening here. Why is it that bgzf_read(fp, x, 4) for example is faster with the bgzf_read_fast inline? Some of it is the removal of a function, but I think a lot of it is because when we inline a memcpy of constant 4 bytes it becomes a straight 32-bit load to register. When we call the function, the length is a variable so it has to call memcpy direct and that in turn has a lot of code in it for dealing with arbitrary length and is probably optimised for copying long things more than short things.

So anywhere we have a bgzf_read with a small-ish constant size, it's generally good practice to call the inline variant so the compiler can optimise away memcpy.

The optimisation can never work in the specific case of bgzf_read_small(fp, x, 32) due to the earlier line:

if (fp->block_length - fp->block_offset > 32) {

i.e. the same optimisation (plus a zero-copy one) happens elsewhere.

Ah yes I see what you mean. It makes no difference to CPU from what I can see (maybe the compiler detects it as an impossibility and removes the code?), but yes it's pointless to go via the static inline.

Now we have the same code for both paths in parsing the x[] array/pointer there's no as much gain to be had on the fp->block_length - fp->block_offset > 32 check anyway (if we kept the "small" variant), however it does still seem to be about 5% faster on small small-read datasets so we may as well keep it.

@16

The bgzf_read function is long and not something that's likely to get inlined even if it was in the header. However most of the time our calls request a small amount of data and they fit within the buffer we've read, so we offer a static inline to do the memcpy when we can, falling back to the long function when we cannot. In terms of CPU time it's not much difference, but the key thing is that it's often CPU time saved in a main thread given the bulk of the decode is often threaded. An example of test_view -B -@16 develop: real 0m48.158s user 6m2.901s sys 0m28.134s real 0m48.730s user 6m3.707s sys 0m28.473s real 0m48.653s user 6m5.215s sys 0m28.637s This PR: real 0m41.731s user 5m59.780s sys 0m30.393s real 0m41.945s user 6m0.367s sys 0m30.426s So we can see it's a consistent win when threading, potentially 10-15% faster throughput depending on work loads.

@16

When running in a high thread count, our the decompression stage in bgzf_read is often the only threaded part meaning whatever is left in main can become the bottleneck once we have sufficient number of threads running. Hence speeding up anything in bam_read1 is key. - sam_realloc_bam_data has an extra 32 bytes. This may not seem much, especially after rounding up to a power of 2. However in tests it makes a significant reduction to memory copies (and also strangely memory size). Tested with both GNU malloc and tcmalloc. - bam_tag2cigar speed up by reordering the checks and simplifying the expression to look for the necessary cigar field. - Avoid a bgzf read and copying from bgzf buffer to the temporary x[] when we know we can copy direct. This subverts the bgzf interface, but it's internal code. Some benchmarks of test_view -B -@16 in.bam develop: 9.73 sec prev commit: 8.32 sec this commit: 7.81 sec Combined this is ~25% speed up.

@32

As we did with bgzf_read_small, shortcutting the big bgzf_write function for the common use case of a small write that fits into the buffer can help reduce the pressure on the main thread. Benchmarks with test_view: Previous commit: for i in `seq 1 3`;do taskset -c 0-31 /usr/bin/time -f '%U user\t%S system\t%e elapsed\t%P %%CPU' ./test/test_view -@32 -b ~/scratch/data/novaseq.10m.bam -p /tmp/_.bam ;done;md5sum /tmp/_.bam 89.81 user 2.22 system 4.11 elapsed 2237% %CPU 89.57 user 2.43 system 4.20 elapsed 2189% %CPU 88.44 user 2.30 system 3.96 elapsed 2291% %CPU bc9ca86ebef3b6669fc7b6fdd7e1acb6 /tmp/_.bam This commit: for i in `seq 1 3`;do taskset -c 0-31 /usr/bin/time -f '%U user\t%S system\t%e elapsed\t%P %%CPU' ./test/test_view -@32 -b ~/scratch/data/novaseq.10m.bam -p /tmp/_.bam ;done;md5sum /tmp/_.bam 86.45 user 1.91 system 3.49 elapsed 2531% %CPU 86.28 user 1.84 system 3.43 elapsed 2562% %CPU 86.81 user 2.19 system 3.54 elapsed 2509% %CPU bc9ca86ebef3b6669fc7b6fdd7e1acb6 /tmp/_.bam So that's about 14% faster throughput. It harms some over places, so this isn't a blanket bgzf_write to bgzf_write_small define. Also following the observation above, similarly restricted bgzf_read_small to bam_read1 instead. The indirection via the small function harms big reads, which affects SAM reading. Benchmarks on ./test/test_view -@32 -b /tmp/_.sam.gz -p /tmp/_.bam shows this speeds up from 3.7s elapsed to 3.4s elapsed. Small, but consistent.

HTSlib starts a new block if an alignment is likely to overflow the current one, so for its own data this will only happen for records longer than 64kbytes. As other implementations may not do this, check that reading works correctly on some BAM files where records have been deliberately split between BGZF blocks. Additionally, check the writing side by making a record with enough CIGAR entries to make it split into multiple BGZF blocks.

daviesrob · 2024-06-04T11:30:32Z

Added a couple of extra tests...

This bug crept in with samtools#1772 which was added since last release, so there is no regression. Fixes samtools#1798 with thanks to John Marshall

Notice: this is the last SAMtools / HTSlib release where CRAM 3.0 will be the default CRAM version. From the next we will change to CRAM 3.1 unless the version is explicitly specified, for example using "samtools view -O cram,version=3.0". Updates ------- * Extend annot-tsv with several new command line options. --delim permits use of other delimiters. --headers for selection of other header formats. --no-header-idx to suppress column index numbers in header. Also removed -h as it is now short for --headers. Note --help still works. (PR samtools#1779) * Allow annot-tsv -a to rename annotations. (PR samtools#1709) * Extend annot-tsv --overlap to be able to specify the overlap fraction separately for source and target. (PR samtools#1811) * Added new APIs to facilitate low-level CRAM container manipulations, used by the new "samtools cat" region filtering code. Functions are: cram_container_get_coords() cram_filter_container() cram_index_extents() cram_container_num2offset() cram_container_offset2num() cram_num_containers() cram_num_containers_between() Also improved cram_index_query() to cope with HTS_IDX_NOCOOR regions. (PR samtools#1771) * Bgzip now retains file modification and access times when compressing and decompressing. (PR samtools#1727, fixes samtools#1718. Requested by Gert Hulselmans.) * Use FNV1a for string hashing in khash. The old algorithm was particularly weak with base-64 style strings and lead to a large number of collisions. (PR samtools#1806. Fixes samtools/samtools#2066, reported by Hans-Joachim Ruscheweyh) * Improve the speed of the nibble2base() function on Intel (PR samtools#1667, PR samtools#1764, PR samtools#1786, PR samtools#1802, thanks to Ruben Vorderman) and ARM (PR samtools#1795, thanks to John Marshall). * bgzf_getline() will now warn if it encounters UTF-16 data. (PR samtools#1487, thanks to John Marshall) * Speed up bgzf_read(). While this does not reduce CPU significantly, it does increase the maximum parallelism available permitting 10-15% faster decoding. (PR samtools#1772, PR samtools#1800, Issue samtools#1798) * Speed up faidx by use of better isgraph methods (PR samtools#1797) and whole-line reading (PR samtools#1799, thanks to John Marshall). * Speed up kputll() function, speeding up BAM -> SAM conversion by about 5% and also samtools depth. (PR samtools#1805) * Added more example code, covering fasta/fastq indexing, tabix indexing and use of the thread pool. (PR samtools#1666) Build Changes ------------- * Code warning fixes for pedantic compilers (PR samtools#1777) and avoid some undefined behaviour (PR samtools#1810, PR samtools#1816, PR samtools#1828). * Windows based CI has been migrated from AppVeyor to GitHub Actions. (PR samtools#1796, PR samtools#1803, PR samtools#1808) * Miscellaneous minor build infrastructure and code fixes. (PR samtools#1807, PR samtools#1829, both thanks to John Marshall) * Updated htscodecs submodule to version 1.6.1 (PR samtools#1828) * Fixed an awk script in the Makefile that only worked with gawk. (PR samtools#1831) Bug fixes --------- * Fix small OSS-Fuzz reported issues with CRAM encoding and long CIGARS and/or illegal positions. (PR samtools#1775, PR samtools#1801, PR samtools#1817) * Fix issues with on-the-fly indexing of VCF/BCF (bcftools --write-index) when not using multiple threads. (PR samtools#1837. Fixes samtools/bcftools#2267, reported by Giulio Genovese) * Stricter limits on POS / MPOS / TLEN in sam_parse1(). This fixes a signed overflow reported by OSS-Fuzz and should help prevent other as-yet undetected bugs. (PR samtools#1812) * Check that the underlying file open worked for preload: URLs. Fixes a NULL pointer dereference reported by OSS-Fuzz. (PR samtools#1821) * Fix an infinite loop in hts_itr_query() when given extremely large positions which cause integer overflow. Also adds hts_bin_maxpos() and hts_idx_maxpos() functions. (PR samtools#1774, thanks to John Marshall and reported by Jesus Alberto Munoz Mesa) * Fix an out of bounds read in hts_itr_multi_next() when switching chromosomes. This bug is present in releases 1.11 to 1.20. (PR samtools#1788. Fixes samtools/samtools#2063, reported by acorvelo) * Work around parsing problems with colons in CHROM names. Fixes samtools/bcftools#2139. (PR samtools#1781, John Marshall / James Bonfield) * Correct the CPU detection for Mac OS X 10.7. cpuid is used by htscodecs (see samtools/htscodecs#116), and the corresponding changes in htslib are PR samtools#1785. Reported by Ryan Carsten Schmidt. * Make BAM zero-length intervals work the same as CRAM; permitted and returning overlapping records. (PR samtools#1787. Fixes samtools/samtools#2060, reported by acorvelo) * Replace assert() with abort() in BCF synced reader. This is not an ideal solution, but it gives consistent behaviour when compiling with or without NDEBUG. (PR samtools#1791, thanks to Martin Pollard) * Fixed failure to change the write block size on compressed SAM or VCF files due to an internal type confusion. (PR samtools#1826) * Fixed an out-of-bounds read in cram_codec_iter_next() (PR samtools#1832)

htslib release 1.21: The primary user-visible changes in this release are updates to the annot-tsv tool and some speed improvements. Full details of other changes and bugs fixed are below. Notice: this is the last SAMtools / HTSlib release where CRAM 3.0 will be the default CRAM version. From the next we will change to CRAM 3.1 unless the version is explicitly specified, for example using "samtools view -O cram,version=3.0". Updates ------- * Extend annot-tsv with several new command line options. --delim permits use of other delimiters. --headers for selection of other header formats. --no-header-idx to suppress column index numbers in header. Also removed -h as it is now short for --headers. Note --help still works. (PR samtools#1779) * Allow annot-tsv -a to rename annotations. (PR samtools#1709) * Extend annot-tsv --overlap to be able to specify the overlap fraction separately for source and target. (PR samtools#1811) * Added new APIs to facilitate low-level CRAM container manipulations, used by the new "samtools cat" region filtering code. Functions are: cram_container_get_coords() cram_filter_container() cram_index_extents() cram_container_num2offset() cram_container_offset2num() cram_num_containers() cram_num_containers_between() Also improved cram_index_query() to cope with HTS_IDX_NOCOOR regions. (PR samtools#1771) * Bgzip now retains file modification and access times when compressing and decompressing. (PR samtools#1727, fixes samtools#1718. Requested by Gert Hulselmans.) * Use FNV1a for string hashing in khash. The old algorithm was particularly weak with base-64 style strings and lead to a large number of collisions. (PR samtools#1806. Fixes samtools/samtools#2066, reported by Hans-Joachim Ruscheweyh) * Improve the speed of the nibble2base() function on Intel (PR samtools#1667, PR samtools#1764, PR samtools#1786, PR samtools#1802, thanks to Ruben Vorderman) and ARM (PR samtools#1795, thanks to John Marshall). * bgzf_getline() will now warn if it encounters UTF-16 data. (PR samtools#1487, thanks to John Marshall) * Speed up bgzf_read(). While this does not reduce CPU significantly, it does increase the maximum parallelism available permitting 10-15% faster decoding. (PR samtools#1772, PR samtools#1800, Issue samtools#1798) * Speed up faidx by use of better isgraph methods (PR samtools#1797) and whole-line reading (PR samtools#1799, thanks to John Marshall). * Speed up kputll() function, speeding up BAM -> SAM conversion by about 5% and also samtools depth. (PR samtools#1805) * Added more example code, covering fasta/fastq indexing, tabix indexing and use of the thread pool. (PR samtools#1666) Build Changes ------------- * Code warning fixes for pedantic compilers (PR samtools#1777) and avoid some undefined behaviour (PR samtools#1810, PR samtools#1816, PR samtools#1828). * Windows based CI has been migrated from AppVeyor to GitHub Actions. (PR samtools#1796, PR samtools#1803, PR samtools#1808) * Miscellaneous minor build infrastructure and code fixes. (PR samtools#1807, PR samtools#1829, both thanks to John Marshall) * Updated htscodecs submodule to version 1.6.1 (PR samtools#1828) * Fixed an awk script in the Makefile that only worked with gawk. (PR samtools#1831) Bug fixes --------- * Fix small OSS-Fuzz reported issues with CRAM encoding and long CIGARS and/or illegal positions. (PR samtools#1775, PR samtools#1801, PR samtools#1817) * Fix issues with on-the-fly indexing of VCF/BCF (bcftools --write-index) when not using multiple threads. (PR samtools#1837. Fixes samtools/bcftools#2267, reported by Giulio Genovese) * Stricter limits on POS / MPOS / TLEN in sam_parse1(). This fixes a signed overflow reported by OSS-Fuzz and should help prevent other as-yet undetected bugs. (PR samtools#1812) * Check that the underlying file open worked for preload: URLs. Fixes a NULL pointer dereference reported by OSS-Fuzz. (PR samtools#1821) * Fix an infinite loop in hts_itr_query() when given extremely large positions which cause integer overflow. Also adds hts_bin_maxpos() and hts_idx_maxpos() functions. (PR samtools#1774, thanks to John Marshall and reported by Jesus Alberto Munoz Mesa) * Fix an out of bounds read in hts_itr_multi_next() when switching chromosomes. This bug is present in releases 1.11 to 1.20. (PR samtools#1788. Fixes samtools/samtools#2063, reported by acorvelo) * Work around parsing problems with colons in CHROM names. Fixes samtools/bcftools#2139. (PR samtools#1781, John Marshall / James Bonfield) * Correct the CPU detection for Mac OS X 10.7. cpuid is used by htscodecs (see samtools/htscodecs#116), and the corresponding changes in htslib are PR samtools#1785. Reported by Ryan Carsten Schmidt. * Make BAM zero-length intervals work the same as CRAM; permitted and returning overlapping records. (PR samtools#1787. Fixes samtools/samtools#2060, reported by acorvelo) * Replace assert() with abort() in BCF synced reader. This is not an ideal solution, but it gives consistent behaviour when compiling with or without NDEBUG. (PR samtools#1791, thanks to Martin Pollard) * Fixed failure to change the write block size on compressed SAM or VCF files due to an internal type confusion. (PR samtools#1826) * Fixed an out-of-bounds read in cram_codec_iter_next() (PR samtools#1832) # -----BEGIN PGP SIGNATURE----- # # iQJJBAABCgAzFiEEaGn1wQ1nqxSLDC2XHsvluhXGTMwFAmbjBlMVHHJtZCtnaXRA # c2FuZ2VyLmFjLnVrAAoJEB7L5boVxkzMgb0P/3AAxkJNLGK8Q1pBBo+in7KzwBOc # /nRfvReki47r4+i1ev98abEF/+XZoYSS8ZkiJH8JMkLCclb1ch2fDQTIL/sjJwrU # MfhYx6oHhzLUOjpG24m390wlvlkEsg5S69jbgImWxLQ1n4ZcSuZGumiwPx5274w5 # 3R/K9vIXtuSF9DTD/jrTgoZeBiN+B4XefAJqHLSo4/XzDo7DNMSZryNGFpHdyKcB # oEfEXX45vJLgzA3Fxo3LqNpC8UTDO5C9u9Uup+hCLN7RTBjV5Fu0UxcwWTSPqfiC # Jroi2FhxNPlDkRL3jraI8LrB3Cqb0G/31OU+cryasrqYnfJrl0e/MyXiRSa0sg25 # uHHWVyyC4JR7uBfUCppTVf76iajbkIVpHgHLZS1heGwM4PCL+AalvsIgGcvLXcBs # g/AwZLoXzao+pbzg8hzN1JuLB5hHcrHp6TsQ6JbTTrWaLeAL3s/+eQ4+d+Np2iUn # v8TMg6onu1HiNgVMQcFXSLC2Q51rB+lD6AA3eusKjcQsaxyk+U7SLCx3OhbvxGsd # SjxQWP1W+WW//e9nm5f4CTSSf1PHoDMyR/kVEc3qwEWJslri68RyB69HrDIVzrWT # 3gKAKJYafi5ssFHsKTsCngESpDsTg3aQV1AGOQPLN/TIjyYDTUOL+/HylRs3k9XN # KOGO5yggx+yqeL3w # =lHVl # -----END PGP SIGNATURE----- # gpg: Signature made Thu 12 Sep 2024 08:18:43 AM PDT # gpg: using RSA key 6869F5C10D67AB148B0C2D971ECBE5BA15C64CCC # gpg: issuer "rmd+git@sanger.ac.uk" # gpg: Can't check signature: No public key

jkbonfield force-pushed the bgzf_inline branch 2 times, most recently from 14d0582 to eca5c35 Compare April 19, 2024 14:41

jkbonfield force-pushed the bgzf_inline branch 2 times, most recently from 27b2c5a to cda6d59 Compare April 29, 2024 13:21

daviesrob self-assigned this May 2, 2024

daviesrob reviewed May 9, 2024

View reviewed changes

jkbonfield force-pushed the bgzf_inline branch from cda6d59 to 5f30c8e Compare May 9, 2024 14:12

jkbonfield added 4 commits May 14, 2024 15:30

Changes from review

84bbe2c

jkbonfield force-pushed the bgzf_inline branch from 5f30c8e to 84bbe2c Compare May 14, 2024 14:30

daviesrob merged commit 11205a9 into samtools:develop Jun 4, 2024

jmarshall mentioned this pull request Jun 28, 2024

Using bgzf_read_small() breaks subsequent bgzf_useek() calls #1798

Closed

jkbonfield mentioned this pull request Jul 1, 2024

Fix a missing bgzf->uncompressed_address incr in bgzf_read_small #1800

Merged

Uh oh!

Conversation

jkbonfield commented Apr 19, 2024

Uh oh!

jkbonfield commented Apr 19, 2024

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jkbonfield May 9, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

daviesrob commented Jun 4, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

jkbonfield May 9, 2024 •

edited

Loading