#define HTS_BUILDING_LIBRARY // Enables HTSLIB_EXPORT, see htslib/hts_defs.h

#include <config.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <errno.h>

#include <unistd.h>

#include <assert.h>

#include <pthread.h>

#include <sys/types.h>

#include <inttypes.h>

#include <zlib.h>

#ifdef HAVE_LIBDEFLATE

#include <libdeflate.h>

#include "htslib/hts.h"

#include "htslib/bgzf.h"

#include "htslib/hfile.h"

#include "htslib/thread_pool.h"

#include "htslib/hts_alloc.h"

#include "htslib/hts_endian.h"

#include "cram/pooled_alloc.h"

#include "hts_internal.h"

#include "bgzf_internal.h"

#include "htslib/khash.h"

#ifndef EFTYPE

#define EFTYPE ENOEXEC

#define BGZF_MT

#define BLOCK_HEADER_LENGTH 18

#define BLOCK_FOOTER_LENGTH 8

static const uint8_t g_magic[19] = "\037\213\010\4\0\0\0\0\0\377\6\0\102\103\2\0\0\0";

typedef struct {

int size;

uint8_t *block;

int64_t end_offset;

} cache_t;

KHASH_MAP_INIT_INT64(bgzf_cache, cache_t)

#ifdef BGZF_MT

typedef struct bgzf_job {

} bgzf_job;

enum mtaux_cmd {

};

typedef struct {

hts_pos_t beg, end;

int tid, is_mapped; // args for hts_idx_push

uint64_t offset, block_number;

} hts_idx_cache_entry;

typedef struct {

int nentries, mentries; // used and allocated

hts_idx_cache_entry *e; // hts_idx elements

} hts_idx_cache_t;

typedef struct bgzf_mtaux_t {

} mtaux_t;

typedef struct

{

uint64_t uaddr; // offset w.r.t. uncompressed data

uint64_t caddr; // offset w.r.t. compressed data

}

bgzidx1_t;

struct bgzidx_t

};

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) {

hts_idx_cache_entry *e;

mtaux_t *mt = fp->mt;

if (!mt)

return hts_idx_push(hidx, tid, beg, end, offset, is_mapped);

// Early check for out of range positions which would fail in hts_idx_push()

if (hts_idx_check_range(hidx, tid, beg, end) < 0)

return -1;

pthread_mutex_lock(&mt->idx_m);

mt->hts_idx = hidx;

hts_idx_cache_t *ic = &mt->idx_cache;

if (ic->nentries >= ic->mentries) {

int new_sz = ic->mentries ? ic->mentries*2 : 1024;

if (!(e = hts_realloc_p(ic->e, sizeof(*ic->e), new_sz))) {

pthread_mutex_unlock(&mt->idx_m);

return -1;

}

ic->e = e;

ic->mentries = new_sz;

}

e = &ic->e[ic->nentries++];

e->tid = tid;

e->beg = beg;

e->end = end;

e->is_mapped = is_mapped;

e->offset = offset & 0xffff;

e->block_number = mt->block_number;

pthread_mutex_unlock(&mt->idx_m);

return 0;

}

static int bgzf_idx_flush(BGZF *fp,

size_t block_uncomp_len, size_t block_comp_len) {

mtaux_t *mt = fp->mt;

if (!mt->idx_cache.e) {

mt->block_written++;

return 0;

}

pthread_mutex_lock(&mt->idx_m);

hts_idx_cache_entry *e = mt->idx_cache.e;

int i;

assert(mt->idx_cache.nentries == 0 || mt->block_written <= e[0].block_number);

for (i = 0; i < mt->idx_cache.nentries && e[i].block_number == mt->block_written; i++) {

if (block_uncomp_len > 0 && e[i].offset == block_uncomp_len) {

/*

// Assert that this is the last entry for the current block_number

assert(i == mt->idx_cache.nentries - 1

|| e[i].block_number < e[i + 1].block_number);

// Work out where the next block starts. For this entry, the

// offset will be zero.

uint64_t next_block_addr = mt->block_address + block_comp_len;

if (hts_idx_push(mt->hts_idx, e[i].tid, e[i].beg, e[i].end,

next_block_addr << 16, e[i].is_mapped) < 0) {

pthread_mutex_unlock(&mt->idx_m);

return -1;

}

// Count this entry and drop out of the loop

i++;

break;

}

if (hts_idx_push(mt->hts_idx, e[i].tid, e[i].beg, e[i].end,

(mt->block_address << 16) + e[i].offset,

e[i].is_mapped) < 0) {

pthread_mutex_unlock(&mt->idx_m);

return -1;

}

}

memmove(&e[0], &e[i], (mt->idx_cache.nentries - i) * sizeof(*e));

mt->idx_cache.nentries -= i;

mt->block_written++;

pthread_mutex_unlock(&mt->idx_m);

return 0;

}

void bgzf_index_destroy(BGZF *fp);

static int bgzf_index_add_block(BGZF *fp);

static int mt_destroy(mtaux_t *mt);

static int bgzf_index_grow(bgzidx_t *idx)

{

if (idx->noffs < idx->moffs)

return 0;

size_t new_moffs = idx->moffs > 0 ? idx->moffs * 2 : 4;

bgzidx1_t *new_offs = hts_realloc_p(idx->offs, sizeof(*new_offs),

new_moffs);

if (new_offs == NULL)

return -1;

idx->offs = new_offs;

idx->moffs = new_moffs;

return 0;

}

static inline void packInt16(uint8_t *buffer, uint16_t value)

{

buffer[0] = value;

buffer[1] = value >> 8;

}

static inline int unpackInt16(const uint8_t *buffer)

{

return buffer[0] | buffer[1] << 8;

}

static inline void packInt32(uint8_t *buffer, uint32_t value)

{

buffer[0] = value;

buffer[1] = value >> 8;

buffer[2] = value >> 16;

buffer[3] = value >> 24;

}

static void razf_info(hFILE *hfp, const char *filename)

{

uint64_t usize, csize;

off_t sizes_pos;

if (filename == NULL || strcmp(filename, "-") == 0) filename = "FILE";

// RAZF files end with USIZE,CSIZE stored as big-endian uint64_t

if ((sizes_pos = hseek(hfp, -16, SEEK_END)) < 0) goto no_sizes;

if (hread(hfp, &usize, 8) != 8 || hread(hfp, &csize, 8) != 8) goto no_sizes;

if (!ed_is_big()) ed_swap_8p(&usize), ed_swap_8p(&csize);

if (csize >= sizes_pos) goto no_sizes; // Very basic validity check

hts_log_error(

" truncate -s %" PRIu64 " %s\n"

"The resulting uncompressed file should be %" PRIu64 " bytes in length.\n"

"likely produce a \"trailing garbage ignored\" message, which can be ignored).",

csize, filename, filename, usize);

return;

no_sizes:

hts_log_error(

"usually be safely ignored.", filename);

}

static const char *bgzf_zerr(int errnum, z_stream *zs)

{

static char buffer[32];

/* Return zs->msg if available.

if (zs && zs->msg) return zs->msg;

// gzerror OF((gzFile file, int *errnum)

switch (errnum) {

case Z_ERRNO:

return strerror(errno);

case Z_STREAM_ERROR:

return "invalid parameter/compression level, or inconsistent stream state";

case Z_DATA_ERROR:

return "invalid or incomplete IO";

case Z_MEM_ERROR:

return "out of memory";

case Z_BUF_ERROR:

return "progress temporarily not possible, or in() / out() returned an error";

case Z_VERSION_ERROR:

return "zlib version mismatch";

case Z_NEED_DICT:

return "data was compressed using a dictionary";

case Z_OK: // 0: maybe gzgets error Z_NULL

default:

snprintf(buffer, sizeof(buffer), "[%d] unknown", errnum);

return buffer; // FIXME: Not thread-safe.

}

}

static BGZF *bgzf_read_init(hFILE *hfpr, const char *filename)

{

BGZF *fp;

uint8_t magic[18];

ssize_t n = hpeek(hfpr, magic, 18);

if (n < 0) return NULL;

fp = (BGZF*)calloc(1, sizeof(BGZF));

if (fp == NULL) return NULL;

fp->is_write = 0;

fp->uncompressed_block = hts_malloc_p(2, BGZF_MAX_BLOCK_SIZE);

if (fp->uncompressed_block == NULL) { free(fp); return NULL; }

fp->compressed_block = (char *)fp->uncompressed_block + BGZF_MAX_BLOCK_SIZE;

fp->is_compressed = (n==18 && magic[0]==0x1f && magic[1]==0x8b);

fp->is_gzip = ( !fp->is_compressed || ((magic[3]&4) && memcmp(&magic[12], "BC\2\0",4)==0) ) ? 0 : 1;

if (fp->is_compressed && (magic[3]&4) && memcmp(&magic[12], "RAZF", 4)==0) {

hts_log_error("Cannot decompress legacy RAZF format");

razf_info(hfpr, filename);

free(fp->uncompressed_block);

free(fp);

errno = EFTYPE;

return NULL;

}

if (!(fp->cache = malloc(sizeof(*fp->cache)))) {

free(fp->uncompressed_block);

free(fp);

return NULL;

}

if (!(fp->cache->h = kh_init(bgzf_cache))) {

free(fp->uncompressed_block);

free(fp->cache);

free(fp);

return NULL;

}

fp->cache->last_pos = 0;

fp->cache->private_data = NULL;

fp->cache->private_data_cleanup = (bgzf_private_data_cleanup_func *) NULL;

return fp;

}

static int mode2level(const char *mode)

{

int i, compress_level = -1;

for (i = 0; mode[i]; ++i)

if (mode[i] >= '0' && mode[i] <= '9') break;

if (mode[i]) compress_level = (int)mode[i] - '0';

if (strchr(mode, 'u')) compress_level = -2;

return compress_level;

}

static BGZF *bgzf_write_init(const char *mode)

{

BGZF *fp;

fp = (BGZF*)calloc(1, sizeof(BGZF));

if (fp == NULL) goto mem_fail;

fp->is_write = 1;

fp->cache = malloc(sizeof(bgzf_cache_t));

if (!fp->cache)

goto mem_fail;

fp->cache->h = NULL;

fp->cache->last_pos = 0;

fp->cache->private_data = NULL;

fp->cache->private_data_cleanup = (bgzf_private_data_cleanup_func *) NULL;

int compress_level = mode2level(mode);

if ( compress_level==-2 )

{

fp->is_compressed = 0;

return fp;

}

fp->is_compressed = 1;

fp->uncompressed_block = hts_malloc_p(2, BGZF_MAX_BLOCK_SIZE);

if (fp->uncompressed_block == NULL) goto mem_fail;

fp->compressed_block = (char *)fp->uncompressed_block + BGZF_MAX_BLOCK_SIZE;

fp->compress_level = compress_level < 0? Z_DEFAULT_COMPRESSION : compress_level; // Z_DEFAULT_COMPRESSION==-1

if (fp->compress_level > 9) fp->compress_level = Z_DEFAULT_COMPRESSION;

if ( strchr(mode,'g') )

{

// gzip output

fp->is_gzip = 1;

fp->gz_stream = (z_stream*)calloc(1,sizeof(z_stream));

if (fp->gz_stream == NULL) goto mem_fail;

fp->gz_stream->zalloc = NULL;

fp->gz_stream->zfree = NULL;

fp->gz_stream->msg = NULL;

int ret = deflateInit2(fp->gz_stream, fp->compress_level, Z_DEFLATED, 15|16, 8, Z_DEFAULT_STRATEGY);

if (ret!=Z_OK) {

hts_log_error("Call to deflateInit2 failed: %s", bgzf_zerr(ret, fp->gz_stream));

goto fail;

}

}

return fp;

mem_fail:

hts_log_error("%s", strerror(errno));

fail:

if (fp != NULL) {

free(fp->cache);

free(fp->uncompressed_block);

free(fp->gz_stream);

free(fp);

}

return NULL;

}

BGZF *bgzf_open(const char *path, const char *mode)

{

BGZF *fp = 0;

if (strchr(mode, 'r')) {

hFILE *fpr;

if ((fpr = hopen(path, mode)) == 0) return 0;

fp = bgzf_read_init(fpr, path);

if (fp == 0) { hclose_abruptly(fpr); return NULL; }

fp->fp = fpr;

} else if (strchr(mode, 'w') || strchr(mode, 'a')) {

hFILE *fpw;

if ((fpw = hopen(path, mode)) == 0) return 0;

fp = bgzf_write_init(mode);

if (fp == NULL) return NULL;

fp->fp = fpw;

}

else { errno = EINVAL; return 0; }

fp->is_be = ed_is_big();

return fp;

}

BGZF *bgzf_dopen(int fd, const char *mode)

{

BGZF *fp = 0;

if (strchr(mode, 'r')) {

hFILE *fpr;

if ((fpr = hdopen(fd, mode)) == 0) return 0;

fp = bgzf_read_init(fpr, NULL);

if (fp == 0) { hclose_abruptly(fpr); return NULL; } // FIXME this closes fd

fp->fp = fpr;

} else if (strchr(mode, 'w') || strchr(mode, 'a')) {

hFILE *fpw;

if ((fpw = hdopen(fd, mode)) == 0) return 0;

fp = bgzf_write_init(mode);

if (fp == NULL) return NULL;

fp->fp = fpw;

}

else { errno = EINVAL; return 0; }

fp->is_be = ed_is_big();

return fp;

}

BGZF *bgzf_hopen(hFILE *hfp, const char *mode)

{

BGZF *fp = NULL;

if (strchr(mode, 'r')) {

fp = bgzf_read_init(hfp, NULL);

if (fp == NULL) return NULL;

} else if (strchr(mode, 'w') || strchr(mode, 'a')) {

fp = bgzf_write_init(mode);

if (fp == NULL) return NULL;

}

else { errno = EINVAL; return 0; }

fp->fp = hfp;

fp->is_be = ed_is_big();

return fp;

}

#ifdef HAVE_LIBDEFLATE

uint32_t hts_crc32(uint32_t crc, const void *buf, size_t len) {

return libdeflate_crc32(crc, buf, len);

}

int bgzf_compress(void *_dst, size_t *dlen, const void *src, size_t slen, int level)

{

if (slen == 0) {

// EOF block

if (*dlen < 28) return -1;

memcpy(_dst, "\037\213\010\4\0\0\0\0\0\377\6\0\102\103\2\0\033\0\3\0\0\0\0\0\0\0\0\0", 28);

*dlen = 28;

return 0;

}

uint8_t *dst = (uint8_t*)_dst;

if (level == 0) {

// Uncompressed data

if (*dlen < slen+5 + BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH) return -1;

dst[BLOCK_HEADER_LENGTH] = 1; // BFINAL=1, BTYPE=00; see RFC1951

u16_to_le(slen, &dst[BLOCK_HEADER_LENGTH+1]); // length

u16_to_le(~slen, &dst[BLOCK_HEADER_LENGTH+3]); // ones-complement length

memcpy(dst + BLOCK_HEADER_LENGTH+5, src, slen);

*dlen = slen+5 + BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH;

} else {

level = level > 0 ? level : 6; // libdeflate doesn't honour -1 as default

// NB levels go up to 12 here.

int lvl_map[] = {0,1,2,3,5,6,7,8,10,12};

level = lvl_map[level>9 ?9 :level];

struct libdeflate_compressor *z = libdeflate_alloc_compressor(level);

if (!z) return -1;

// Raw deflate

size_t clen =

libdeflate_deflate_compress(z, src, slen,

dst + BLOCK_HEADER_LENGTH,

*dlen - BLOCK_HEADER_LENGTH - BLOCK_FOOTER_LENGTH);

if (clen <= 0) {

hts_log_error("Call to libdeflate_deflate_compress failed");

libdeflate_free_compressor(z);

return -1;

}

*dlen = clen + BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH;

libdeflate_free_compressor(z);

}

// write the header

memcpy(dst, g_magic, BLOCK_HEADER_LENGTH); // the last two bytes are a place holder for the length of the block

packInt16(&dst[16], *dlen - 1); // write the compressed length; -1 to fit 2 bytes

// write the footer

uint32_t crc = libdeflate_crc32(0, src, slen);

packInt32((uint8_t*)&dst[*dlen - 8], crc);

packInt32((uint8_t*)&dst[*dlen - 4], slen);

return 0;

}

uint32_t hts_crc32(uint32_t crc, const void *buf, size_t len) {

return crc32(crc, buf, len);

}

int bgzf_compress(void *_dst, size_t *dlen, const void *src, size_t slen, int level)

{

uint32_t crc;

z_stream zs;

uint8_t *dst = (uint8_t*)_dst;

if (level == 0) {

uncomp:

// Uncompressed data

if (*dlen < slen+5 + BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH) return -1;

dst[BLOCK_HEADER_LENGTH] = 1; // BFINAL=1, BTYPE=00; see RFC1951

u16_to_le(slen, &dst[BLOCK_HEADER_LENGTH+1]); // length

u16_to_le(~slen, &dst[BLOCK_HEADER_LENGTH+3]); // ones-complement length

memcpy(dst + BLOCK_HEADER_LENGTH+5, src, slen);

*dlen = slen+5 + BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH;

} else {

// compress the body

zs.zalloc = NULL; zs.zfree = NULL;

zs.msg = NULL;

zs.next_in = (Bytef*)src;

zs.avail_in = slen;

zs.next_out = dst + BLOCK_HEADER_LENGTH;

zs.avail_out = *dlen - BLOCK_HEADER_LENGTH - BLOCK_FOOTER_LENGTH;

int ret = deflateInit2(&zs, level, Z_DEFLATED, -15, 8, Z_DEFAULT_STRATEGY); // -15 to disable zlib header/footer

if (ret!=Z_OK) {

hts_log_error("Call to deflateInit2 failed: %s", bgzf_zerr(ret, &zs));

return -1;

}

if ((ret = deflate(&zs, Z_FINISH)) != Z_STREAM_END) {

if (ret == Z_OK && zs.avail_out == 0) {

deflateEnd(&zs);

goto uncomp;

} else {

hts_log_error("Deflate operation failed: %s", bgzf_zerr(ret, ret == Z_DATA_ERROR ? &zs : NULL));

}

return -1;

}

// If we used up the entire output buffer, then we either ran out of

// room or we *just* fitted, but either way we may as well store

// uncompressed for faster decode.

if (zs.avail_out == 0) {

deflateEnd(&zs);

goto uncomp;

}

if ((ret = deflateEnd(&zs)) != Z_OK) {

hts_log_error("Call to deflateEnd failed: %s", bgzf_zerr(ret, NULL));

return -1;

}

*dlen = zs.total_out + BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH;

}

// write the header

memcpy(dst, g_magic, BLOCK_HEADER_LENGTH); // the last two bytes are a place holder for the length of the block

packInt16(&dst[16], *dlen - 1); // write the compressed length; -1 to fit 2 bytes

// write the footer

crc = crc32(crc32(0L, NULL, 0L), (Bytef*)src, slen);

packInt32((uint8_t*)&dst[*dlen - 8], crc);

packInt32((uint8_t*)&dst[*dlen - 4], slen);

return 0;

}

#endif // HAVE_LIBDEFLATE

static int bgzf_gzip_compress(BGZF *fp, void *_dst, size_t *dlen, const void *src, size_t slen, int level)

{

uint8_t *dst = (uint8_t*)_dst;

z_stream *zs = fp->gz_stream;

int flush = slen ? Z_PARTIAL_FLUSH : Z_FINISH;

zs->next_in = (Bytef*)src;

zs->avail_in = slen;

zs->next_out = dst;

zs->avail_out = *dlen;

int ret = deflate(zs, flush);

if (ret == Z_STREAM_ERROR) {

hts_log_error("Deflate operation failed: %s", bgzf_zerr(ret, NULL));

return -1;

}

if (zs->avail_in != 0) {

hts_log_error("Deflate block too large for output buffer");

return -1;

}

*dlen = *dlen - zs->avail_out;

return 0;

}

static int deflate_block(BGZF *fp, int block_length)

{

size_t comp_size = BGZF_MAX_BLOCK_SIZE;

int ret;

if ( !fp->is_gzip )

ret = bgzf_compress(fp->compressed_block, &comp_size, fp->uncompressed_block, block_length, fp->compress_level);

else

ret = bgzf_gzip_compress(fp, fp->compressed_block, &comp_size, fp->uncompressed_block, block_length, fp->compress_level);

if ( ret != 0 )

{

hts_log_debug("Compression error %d", ret);

fp->errcode |= BGZF_ERR_ZLIB;

return -1;

}

fp->block_offset = 0;

return comp_size;

}

#ifdef HAVE_LIBDEFLATE

static int bgzf_uncompress(uint8_t *dst, size_t *dlen,

const uint8_t *src, size_t slen,

uint32_t expected_crc) {

struct libdeflate_decompressor *z = libdeflate_alloc_decompressor();

if (!z) {

hts_log_error("Call to libdeflate_alloc_decompressor failed");

return -1;

}

int ret = libdeflate_deflate_decompress(z, src, slen, dst, *dlen, dlen);

libdeflate_free_decompressor(z);

if (ret != LIBDEFLATE_SUCCESS) {

hts_log_error("Inflate operation failed: %d", ret);

return -1;

}

uint32_t crc = libdeflate_crc32(0, (unsigned char *)dst, *dlen);

#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION

// Pretend the CRC was OK so the fuzzer doesn't have to get it right

crc = expected_crc;

if (crc != expected_crc) {

hts_log_error("CRC32 checksum mismatch");

return -2;

}

return 0;

}

static int bgzf_uncompress(uint8_t *dst, size_t *dlen,

const uint8_t *src, size_t slen,

uint32_t expected_crc) {

z_stream zs = {

.zalloc = NULL,

.zfree = NULL,

.msg = NULL,

.next_in = (Bytef*)src,

.avail_in = slen,

.next_out = (Bytef*)dst,

.avail_out = *dlen

};

int ret = inflateInit2(&zs, -15);

if (ret != Z_OK) {

hts_log_error("Call to inflateInit2 failed: %s", bgzf_zerr(ret, &zs));

return -1;

}

if ((ret = inflate(&zs, Z_FINISH)) != Z_STREAM_END) {

hts_log_error("Inflate operation failed: %s", bgzf_zerr(ret, ret == Z_DATA_ERROR ? &zs : NULL));

if ((ret = inflateEnd(&zs)) != Z_OK) {

hts_log_warning("Call to inflateEnd failed: %s", bgzf_zerr(ret, NULL));

}

return -1;

}

if ((ret = inflateEnd(&zs)) != Z_OK) {

hts_log_error("Call to inflateEnd failed: %s", bgzf_zerr(ret, NULL));

return -1;

}

*dlen = *dlen - zs.avail_out;

uint32_t crc = crc32(crc32(0L, NULL, 0L), (unsigned char *)dst, *dlen);

#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION

// Pretend the CRC was OK so the fuzzer doesn't have to get it right

crc = expected_crc;

if (crc != expected_crc) {

hts_log_error("CRC32 checksum mismatch");

return -2;

}

return 0;

}

#endif // HAVE_LIBDEFLATE

static int inflate_block(BGZF* fp, int block_length)

{

size_t dlen = BGZF_MAX_BLOCK_SIZE;

uint32_t crc = le_to_u32((uint8_t *)fp->compressed_block + block_length-8);

int ret = bgzf_uncompress(fp->uncompressed_block, &dlen,

(Bytef*)fp->compressed_block + 18,

block_length - 18, crc);

if (ret < 0) {

if (ret == -2)

fp->errcode |= BGZF_ERR_CRC;

else

fp->errcode |= BGZF_ERR_ZLIB;

return -1;

}

return dlen;

}

static int inflate_gzip_block(BGZF *fp)

{

// we will set this to true when we detect EOF, so we don't bang against the EOF more than once per call

int input_eof = 0;

// write to the part of the output buffer after block_offset

fp->gz_stream->next_out = (Bytef*)fp->uncompressed_block + fp->block_offset;

fp->gz_stream->avail_out = BGZF_MAX_BLOCK_SIZE - fp->block_offset;

while ( fp->gz_stream->avail_out != 0 ) {

// until we fill the output buffer (or hit EOF)

if ( !input_eof && fp->gz_stream->avail_in == 0 ) {

// we are out of input data in the buffer. Get more.

fp->gz_stream->next_in = fp->compressed_block;

int ret = hread(fp->fp, fp->compressed_block, BGZF_BLOCK_SIZE);

if ( ret < 0 ) {

// hread had an error. Pass it on.

return ret;

}

fp->gz_stream->avail_in = ret;

if ( fp->gz_stream->avail_in < BGZF_BLOCK_SIZE ) {

// we have reached EOF but the decompressor hasn't necessarily

input_eof = 1;

}

}

fp->gz_stream->msg = NULL;

// decompress as much data as we can

int ret = inflate(fp->gz_stream, Z_SYNC_FLUSH);

if ( (ret < 0 && ret != Z_BUF_ERROR) || ret == Z_NEED_DICT ) {

// an error occurred, other than running out of space

hts_log_error("Inflate operation failed: %s", bgzf_zerr(ret, ret == Z_DATA_ERROR ? fp->gz_stream : NULL));

fp->errcode |= BGZF_ERR_ZLIB;

return -1;

} else if ( ret == Z_STREAM_END ) {

// we finished a GZIP member

// scratch for peeking to see if the file is over

char c;

if (fp->gz_stream->avail_in > 0 || hpeek(fp->fp, &c, 1) == 1) {

// there is more data; try and read another GZIP member in the remaining data

int reset_ret = inflateReset(fp->gz_stream);

if (reset_ret != Z_OK) {

hts_log_error("Call to inflateReset failed: %s", bgzf_zerr(reset_ret, NULL));

fp->errcode |= BGZF_ERR_ZLIB;

return -1;

}

} else {

// we consumed all the input data and hit Z_STREAM_END

// so stop looping, even if we never fill the output buffer

break;

}

} else if ( ret == Z_BUF_ERROR && input_eof && fp->gz_stream->avail_out > 0 ) {

// the gzip file has ended prematurely

hts_log_error("Gzip file truncated");

fp->errcode |= BGZF_ERR_IO;

return -1;

}

}

// when we get here, the buffer is full or there is an EOF after a complete gzip member

return BGZF_MAX_BLOCK_SIZE - fp->gz_stream->avail_out;

}

static int check_header(const uint8_t *header)

{

if ( header[0] != 31 || header[1] != 139 || header[2] != 8 ) return -2;

return ((header[3] & 4) != 0

&& unpackInt16((uint8_t*)&header[10]) == 6

&& header[12] == 'B' && header[13] == 'C'

&& unpackInt16((uint8_t*)&header[14]) == 2) ? 0 : -1;

}

static void free_cache(BGZF *fp)

{

khint_t k;

if (fp->cache->h) {

khash_t(bgzf_cache) *h = fp->cache->h;

for (k = kh_begin(h); k < kh_end(h); ++k)

if (kh_exist(h, k)) free(kh_val(h, k).block);

kh_destroy(bgzf_cache, h);

}

bgzf_clear_private_data(fp);

free(fp->cache);

}

static int load_block_from_cache(BGZF *fp, int64_t block_address)

{

khint_t k;

cache_t *p;

khash_t(bgzf_cache) *h = fp->cache->h;

k = kh_get(bgzf_cache, h, block_address);

if (k == kh_end(h)) return 0;

p = &kh_val(h, k);

if (fp->block_length != 0) fp->block_offset = 0;

fp->block_address = block_address;

fp->block_length = p->size;

memcpy(fp->uncompressed_block, p->block, p->size);

if ( hseek(fp->fp, p->end_offset, SEEK_SET) < 0 )

{

// todo: move the error up

hts_log_error("Could not hseek to %" PRId64, p->end_offset);

exit(1);

}

return p->size;

}

static void cache_block(BGZF *fp, int size)

{

int ret;

khint_t k, k_orig;

uint8_t *block = NULL;

cache_t *p;

//fprintf(stderr, "Cache block at %llx\n", (int)fp->block_address);

khash_t(bgzf_cache) *h = fp->cache->h;

if (BGZF_MAX_BLOCK_SIZE >= fp->cache_size) return;

if (fp->block_length < 0 || fp->block_length > BGZF_MAX_BLOCK_SIZE) return;

if ((kh_size(h) + 1) * BGZF_MAX_BLOCK_SIZE > (uint32_t)fp->cache_size) {

/* Remove uniformly from any position in the hash by a simple

if (fp->cache->last_pos >= kh_end(h)) fp->cache->last_pos = kh_begin(h);

k_orig = k = fp->cache->last_pos;

if (++k >= kh_end(h)) k = kh_begin(h);

while (k != k_orig) {

if (kh_exist(h, k))

break;

if (++k == kh_end(h))

k = kh_begin(h);

}

fp->cache->last_pos = k;

if (k != k_orig) {

block = kh_val(h, k).block;

kh_del(bgzf_cache, h, k);

}

} else {

block = (uint8_t*)malloc(BGZF_MAX_BLOCK_SIZE);

}

if (!block) return;

k = kh_put(bgzf_cache, h, fp->block_address, &ret);

if (ret <= 0) { // kh_put failed, or in there already (shouldn't happen)

free(block);

return;

}

p = &kh_val(h, k);

p->size = fp->block_length;

p->end_offset = fp->block_address + size;

p->block = block;