Htscodecs
=========

See the NEWS file for a list of updates and version details.

[![Build Status](https://api.cirrus-ci.com/github/jkbonfield/htscodecs.svg?branch=master)](https://cirrus-ci.com/github/jkbonfield/htscodecs)

This repository implements the custom CRAM codecs used for "EXTERNAL"
block types.  These consist of two variants of the rANS codec (8-bit
and 16-bit renormalisation, with run-length encoding and bit-packing
also supported in the latter), a dynamic arithmetic coder, and custom
codecs for name/ID compression and quality score compression derived
from fqzcomp.


They come with small command line test tools to act as both
compression exploration programs and as part of the test harness.


Building
--------

If building from git, you'll need to recreate the configure script
using autoconf.  "autoreconf -i" should work if you have the
appropriate tools.

From then on, it follows the normal "./configure; make" or
"mkdir build; cd build; ../configure; make" rule.

The library can be used as a git sub-module or as a completely
separate entity.  If you are attempting to make use of these codecs
within your own library, such as we do within Staden io_lib, it may be
useful to configure this with `--disable-shared --with-pic'.

Intel GDS/Downfall
------------------

The new Intel microcode to patch the Gather Data Sampling (GDS) bug
has severely impacted the performance of the AVX2 and AVX512 gather
SIMD intrinsics.  So much so that simulating vector gather
instructions using traditional scalar code is now more performant.
Due to the slower AVX512 gathers on AMD Zen4, the simulated gather
code is faster there too, even without having a Downfall mitigation.

Unfortunately this simulated code is slower than the original gather
code when run on an unpatched CPUs.  Without the GDS microcode patch
on Intel CPUs, the simulated gathers slows the rate of decode by
10-30% and the rate of encode for AVX512 by up to 50%.  The GDS
microcode patch however can slow AVX2 decode and AVX512 encode by
2-3x.

Ideally we would auto-detect and select the optimal algorithm, but the
easy solution is simply to pick the method with the most uniform
performance without having bad worst-case rates.

Hence by default the simulated gather implementation is now the
default.  The old code can be reenabled by building with:

    make CPPFLAGS=-DUSE_GATHER


Testing
-------

There is a "make check" rule.  If you're using a modern clang you can
also cd to the tests directory and do "make fuzz" to build some fuzz
testing targets, but you'll likely need to modify Makefile.am first as
this has some hard-coded local paths.

We also provide test data and some command line tools to demonstrate
usage of the compression codecs.  These are in the tests directory
also. Example usage:

    ./fqzcomp_qual -s 1 < dat/q40+dir > /tmp/q40.comp
    ./fqzcomp_qual  -d < /tmp/q40.comp > /tmp/q40.uncomp
    awk '{print $1}' dat/q40+dir | md5sum;    # f91473032dd6942e72abec0868f17161
    awk '{print $1}' /tmp/q40.uncomp | md5sum;# f91473032dd6942e72abec0868f17161

The fqzcomp test format is one quality values per line, with an
optional additional parameter (0 or 1) to indicate READ1 or READ2 flag
status.

There is a larger set of test data in the htscodecs-corpus repository
(https://github.com/jkbonfield/htscodecs-corpus).  If this is cloned
into the tests subdirectory of htscodecs then the htscodecs "make
check" will also use that larger data set for testing purposes.


API
---

Many functions just take an input buffer and size and return an output
buffer, setting *out_size with the decoded size.  NULL is returned for
error.  This buffer is malloced and is expected to be freed by the
caller.  These are the *`compress` and *`uncompress` functions.

A second variant sometimes exists where the output buffer is
optionally allocated by the caller (it may be NULL in which case it
has the same operation as above).  If specified, `*out_size` must also
be set to the allocated size of `out`.  These are the `compress_to`
and `uncompress_to` functions.

The compress size sometimes needs additional options.  For the rANS
and arithmetic coder this is the "order".  Values of 0 and 1 are
simple order-0 and order-1 entropy encoder, but this is a bit field
and the more advanced codecs have additional options to pass in order
(so it should really be renamed to flags).  See below.  Fqzcomp
requires more input data - also see below.  In all cases, sufficient
information is stored in the compressed byte stream such that the
decompression will work without needing these input paramaters.

Finally the various `compress_bound` functions give the size of buffer
needed to be allocated when compressing a block of data.


### Static rANS 4x8 (introduced in CRAM v3.0)

```
#include "htscodecs/rANS_static.h"

unsigned char *rans_compress(unsigned char *in, unsigned int in_size,
                             unsigned int *out_size, int order);
unsigned char *rans_uncompress(unsigned char *in, unsigned int in_size,
                               unsigned int *out_size);
```

This is the earlier rANS entropy encoder using 4 rANS states and 8-bit
renormalisation, with Order-0 and Order-1 entropy models.

No (un)compress_to functions exist for this older codec.


### Static rANS 4x16 and 32x16 with bit-pack/RLE (CRAM v3.1):

```
#include "htscodecs/rANS_static4x16.h"

#define RANS_ORDER_X32    0x04  // 32-way unrolling instead of 4-way
#define RANS_ORDER_STRIPE 0x08  // N streams for every Nth byte (N==order>>8)
#define RANS_ORDER_NOSZ   0x10  // Don't store the original size
#define RANS_ORDER_CAT    0x20  // Nop; for tiny data segments
#define RANS_ORDER_RLE    0x40  // Run length encoding
#define RANS_ORDER_PACK   0x80  // Pack 2,4,8 or infinite symbols into a byte.

unsigned int rans_compress_bound_4x16(unsigned int size, int order);
unsigned char *rans_compress_to_4x16(unsigned char *in,  unsigned int in_size,
                                     unsigned char *out, unsigned int *out_size,
                                     int order);
unsigned char *rans_compress_4x16(unsigned char *in, unsigned int in_size,
                                  unsigned int *out_size, int order);
unsigned char *rans_uncompress_to_4x16(unsigned char *in,  unsigned int in_size,
                                       unsigned char *out, unsigned int *out_size);
unsigned char *rans_uncompress_4x16(unsigned char *in, unsigned int in_size,
                                    unsigned int *out_size);
```

This is a faster version with 16-bit renormalisation and optional
transforms (RLE, small alphabet bit-packing, and interleaving of N
streams for e.g. 32-bit integer compression).  Additionally the
`order` field may include bit `RANS_ORDER_X32` in which case a 32-way
unrolled version will be used instead, with automatic CPU detection
and dispatching to an appropriate SIMD implementation if available.

### Adaptive arithmetic coding (CRAM v3.1):

```
#include "htscodecs/arith_dynamic.h"

unsigned char *arith_compress(unsigned char *in, unsigned int in_size,
                              unsigned int *out_size, int order);

unsigned char *arith_uncompress(unsigned char *in, unsigned int in_size,
                                unsigned int *out_size);

unsigned char *arith_compress_to(unsigned char *in,  unsigned int in_size,
                                 unsigned char *out, unsigned int *out_size,
                                 int order);

unsigned char *arith_uncompress_to(unsigned char *in, unsigned int in_size,
                                   unsigned char *out, unsigned int *out_sz);

unsigned int arith_compress_bound(unsigned int size, int order);
```

These reuse the same `RANS_ORDER` bit fields and abilities above with
the exception of X32 as there is currently no unrolling of this code.

### Name tokeniser (CRAM v3.1):

```
#include "htscodecs/tokenise_name3.h"

uint8_t *encode_names(char *blk, int len, int level, int use_arith,
                      int *out_len, int *last_start_p);

uint8_t *decode_names(uint8_t *in, uint32_t sz, uint32_t *out_len);
```

This differs to the general purpose entropy encoders as it takes a
specific type of data.  The names should be newline or nul separated
for `encode_names`.  `decode_names` will alway return nul terminated
names, so you may need to swap these to newlines if you do round-trip
tests.

The compression level controls how hard it tries to find the optimum
compression method per internal token column.  By default it'll use
the rANS 4x16 codec, but with non-zero `use_arith` it'll use the
adaptive arithmetic coder instead.

If non-NULL, last_start_p can be used to point to a partial name if an
arbitrary block of names were supplied that don't end of a whole read
name. (Is this useful?  Probably not.)


### FQZComp Qual (CRAM v3.1):


```
#include "htscodecs/fqzcomp_qual.h"

#define FQZ_FREVERSE 16
#define FQZ_FREAD2 128

typedef struct {
    int num_records;
    uint32_t *len;    // of size num_records
    uint32_t *flags;  // of size num_records
} fqz_slice;

char *fqz_compress(int vers, fqz_slice *s, char *in, size_t uncomp_size,
                   size_t *comp_size, int strat, fqz_gparams *gp);
char *fqz_decompress(char *in, size_t comp_size, size_t *uncomp_size,
                     int *lengths, int nlengths);
```

This is derived from the quality compression in fqzcomp.  The input
buffer is a concatenated block of quality strings, without any
separator.  In order to achieve maximum compression it needs to know
where these separators are, so they must be passed in via the
`fqz_slice` struct.

The summation of length fields should match the input uncomp_size
field.  Note the len fields may not actually be the length of the
original sequences as some CRAM features may additional quality values
(eg the "B" feature).

It can also be beneficial to supply per-record flags so fqzcomp can
determine whether orientation (complement strand) helps and whether
the READ1 vs READ2 quality distributions differ.  These are just
sub-fields from BAM FLAG.

The fqz_gparams will normally be passed in as NULL and the encoder
will automatically select parameters.  If you wish to fine tune the
compression methods, see the fqz_params and fqz_gparams structures in
the header file.  You may also find the fqz_qual_stats() utility
function helpful for gathering statistics on your quality values.

For decompression, the lengths array is optional and may be specified
as NULL.  If passed in, it must be of size nlengths and it will be
filled out with the decoded length of each quality string.  Note
regardless of whether lengths is NULL or not, the buffer returned will
be concatenated values so there is no way to tell where one record
finishes and the next starts.  (CRAM itself knows this via other means.)