diff --git a/src/slow5_press.c b/src/slow5_press.c index 4997018..d970dc6 100644 --- a/src/slow5_press.c +++ b/src/slow5_press.c @@ -1253,483 +1253,17 @@ static uint32_t *delta_increasing_u32(const uint32_t *in, uint64_t nin) return out; } -static uint32_t get_max_u32(const uint32_t *in, uint64_t nin) -{ - uint32_t max; - uint64_t i; - - max = 0; - - for (i = 0; i < nin; i++) { - if (in[i] > max) - max = in[i]; - } - - return max; -} - - -#define SLOW5_DIV_ROUND_UP(n, d) ((n - 1) / d + 1) -#define SLOW5_BITS_PER_BYTE (8) -#define SLOW5_BITS_TO_BYTES(n) SLOW5_DIV_ROUND_UP(n, SLOW5_BITS_PER_BYTE) -#define SLOW5_BYTES_TO_BITS(n) ((n) * SLOW5_BITS_PER_BYTE) - -/* uintx */ - -static uint64_t uintx_bound(uint8_t in_bits, uint8_t out_bits, uint64_t nin) -{ - uint64_t nin_elems; - - nin_elems = SLOW5_BYTES_TO_BITS(nin) / in_bits; - - if (!out_bits) - return 0; - - return SLOW5_BITS_TO_BYTES(nin_elems * out_bits); -} - - -static void uintx_update(uint8_t in_bits, uint8_t out_bits, uint64_t *in_i, - uint64_t *out_i, uint8_t *in_bits_free, - uint8_t *out_bits_free, uint8_t *bits_left) -{ - int8_t gap; - - gap = in_bits - out_bits; - - if (gap > 0) { - if (gap > *in_bits_free) { - (*in_i)++; - gap -= *in_bits_free; - } - *bits_left = out_bits; - *in_bits_free = SLOW5_BITS_PER_BYTE - gap % SLOW5_BITS_PER_BYTE; - *in_i += gap / SLOW5_BITS_PER_BYTE; - } else { - gap *= -1; - if (gap > *out_bits_free) { - (*out_i)++; - gap -= *out_bits_free; - } - *bits_left = in_bits; - *out_bits_free = SLOW5_BITS_PER_BYTE - gap % SLOW5_BITS_PER_BYTE; - *out_i += gap / SLOW5_BITS_PER_BYTE; - } -} - -/* if in_bits > out_bits: in must be in big endian format */ -static int uintx_press_core(uint8_t in_bits, uint8_t out_bits, const uint8_t *in, - uint64_t nin, uint8_t *out, uint64_t *nout) -{ - /* - * in_bits = 16 - * out_bits = 11 - * gap = 5 - * in_free_bits = 8 - * in = [00000{010, 10001011}, ...] - * out = [{01010001, 011}..., ...] - * - * press operations on P11: - * out[0] = in[0] << 5 | in[1] >> 3; - * out[1] = in[1] << 5 | in[2] >> 3; - * ... - * - * in_bits = 14 - * out_bits = 9 - * gap = 5 - * in_free_bits = 8,2 - * in = [00000{010, 100010}00, 000 ...] - * out = [{01010001, 0}..., ...] - * - * reverse: - * in_bits = 11 - * out_bits = 16 - * gap = -5 - * in = [{01010001, 011}..., ...] - * out = [00000{010, 10001011}, ...] - * - * press operations on compressed P11: - * out[0] = in[0] >> 5 - * out[1] = in[0] << 3 | in[1] >> 5; - * out[2] = in[1] << 3 | in[2] >> 5; - * ... - */ - - int dirty; - int8_t gap; - uint64_t i; - uint64_t in_i; - uint64_t out_i; - uint8_t bits_left; - uint8_t cur_out; - uint8_t in_bits_free; - uint8_t mask; - uint8_t out_bits_free; - - /* when decompressing some of out may be skipped */ - if (in_bits < out_bits) - (void) memset(out, 0, *nout); - - bits_left = 0; - cur_out = 0; - dirty = 0; - i = 0; - in_bits_free = SLOW5_BITS_PER_BYTE; - in_i = 0; - out_bits_free = SLOW5_BITS_PER_BYTE; - out_i = 0; - - uintx_update(in_bits, out_bits, &in_i, &out_i, &in_bits_free, - &out_bits_free, &bits_left); - - while (i < nin) { - mask = 0xFF >> (SLOW5_BITS_PER_BYTE - in_bits_free); - gap = in_bits_free - out_bits_free; - if (gap > 0) { - cur_out |= (in[in_i] & mask) >> gap; - in_bits_free -= out_bits_free; - bits_left -= out_bits_free; - out_bits_free = 0; - } else { - cur_out |= (in[in_i] & mask) << (-1 * gap); - out_bits_free -= in_bits_free; - bits_left -= in_bits_free; - in_bits_free = 0; - } - dirty = 1; - - if (!in_bits_free) { - in_bits_free = SLOW5_BITS_PER_BYTE; - in_i++; - } - - if (!out_bits_free) { - if (out_i == *nout) - return -1; - out[out_i++] = cur_out; - out_bits_free = SLOW5_BITS_PER_BYTE; - cur_out = 0; - dirty = 0; - } - - if (!bits_left) { - uintx_update(in_bits, out_bits, &in_i, &out_i, - &in_bits_free, &out_bits_free, - &bits_left); - i++; - } - } - - /* if there is still data to flush */ - if (dirty) { - if (out_i == *nout) - return -1; - out[out_i++] = cur_out; - } - - *nout = out_i; - - return 0; -} - -static inline uint16_t uint16_bswap(uint16_t x) { - return (x<<8) | (x>>8); -} - -static inline uint32_t uint32_bswap(uint32_t x) { - return (x>>24) | (x>>8 & 0xff00) | (x<<8 & 0xff0000) | (x<<24); -} - -static inline uint64_t uint64_bswap(uint64_t x){ - return ((uint32_bswap(x)+0ULL) << 32) | uint32_bswap(x>>32); -} - -#define slow5_letobe16(x) uint16_bswap(x) -#define slow5_be16tole(x) uint16_bswap(x) -#define slow5_letobe32(x) uint32_bswap(x) -#define slow5_be32tole(x) uint32_bswap(x) -#define slow5_letobe64(x) uint64_bswap(x) -#define slow5_be64tole(x) uint64_bswap(x) - -/* copy h from host endian into be as big endian */ -static void uintx_htobe(uint8_t in_bits, const uint8_t *h, uint8_t *be, uint64_t n) -{ - uint64_t i; - - if (in_bits <= SLOW5_BYTES_TO_BITS(sizeof (uint8_t))) { - return; - } else if (in_bits <= SLOW5_BYTES_TO_BITS(sizeof (uint16_t))) { - for (i = 0; i < n / sizeof (uint16_t); i++) { - ((uint16_t *) be)[i] = slow5_letobe16(((uint16_t *) h)[i]); - } - } else if (in_bits <= SLOW5_BYTES_TO_BITS(sizeof (uint32_t))) { - for (i = 0; i < n / sizeof (uint32_t); i++) { - ((uint32_t *) be)[i] = slow5_letobe32(((uint32_t *) h)[i]); - } - } else if (in_bits <= SLOW5_BYTES_TO_BITS(sizeof (uint64_t))) { - for (i = 0; i < n / sizeof (uint64_t); i++) { - ((uint64_t *) be)[i] = slow5_letobe64(((uint64_t *) h)[i]); - } - } -} - -/* copy be from big endian into h as host endian */ -static void uintx_betoh(uint8_t out_bits, const uint8_t *be, uint8_t *h, uint64_t n) -{ - uint64_t i; - - if (out_bits <= SLOW5_BYTES_TO_BITS(sizeof (uint8_t))) { - return; - } else if (out_bits <= SLOW5_BYTES_TO_BITS(sizeof (uint16_t))) { - for (i = 0; i < n / sizeof (uint16_t); i++) { - ((uint16_t *) h)[i] = slow5_be16tole(((uint16_t *) be)[i]); - } - } else if (out_bits <= SLOW5_BYTES_TO_BITS(sizeof (uint32_t))) { - for (i = 0; i < n / sizeof (uint32_t); i++) { - ((uint32_t *) h)[i] = slow5_be32tole(((uint32_t *) be)[i]); - } - } else if (out_bits <= SLOW5_BYTES_TO_BITS(sizeof (uint64_t))) { - for (i = 0; i < n / sizeof (uint64_t); i++) { - ((uint64_t *) h)[i] = slow5_be64tole(((uint64_t *) be)[i]); - } - } -} - - -static int uintx_press(uint8_t in_bits, uint8_t out_bits, const uint8_t *in, - uint64_t nin, uint8_t *out, uint64_t *nout) -{ - int ret; - uint64_t nin_bytes; - uint8_t *in_be; - - if (!out_bits) { - *nout = 0; - return 0; - } - - nin_bytes = nin * SLOW5_BITS_TO_BYTES(in_bits); - in_be = malloc(nin_bytes); - if (!in_be) - return -1; - uintx_htobe(in_bits, in, in_be, nin_bytes); - - ret = uintx_press_core(in_bits, out_bits, in_be, nin, out, nout); - free(in_be); - - return ret; -} - - -static int uint0_depress(uint8_t out_bits, uint64_t nin, uint8_t *out, uint64_t *nout) -{ - uint64_t nout_bytes; - - nout_bytes = nin * SLOW5_BITS_TO_BYTES(out_bits); - if (*nout < nout_bytes) - return -1; - - (void) memset(out, 0, nout_bytes); - *nout = nout_bytes; - - return 0; -} - - - -static int uintx_depress(uint8_t in_bits, uint8_t out_bits, const uint8_t *in, - uint64_t nin, uint8_t *out, uint64_t *nout) -{ - int ret; - uint8_t *out_be; - - if (!in_bits) - return uint0_depress(out_bits, nin, out, nout); - - out_be = malloc(*nout); - if (!out_be) - return -1; - - ret = uintx_press_core(in_bits, out_bits, in, nin, out_be, nout); - - if (ret == 0) - uintx_betoh(out_bits, out_be, out, *nout); - - free(out_be); - return ret; -} - -#define SLOW5_DEFINE_UINTX(bits) \ -static uint64_t uintx_bound_##bits(uint8_t out_bits, uint64_t nin) \ -{ \ - return uintx_bound(bits, out_bits, nin); \ -} \ -static int uintx_press_##bits(uint8_t out_bits, const uint8_t *in, uint64_t nin, \ - uint8_t *out, uint64_t *nout) \ -{ \ - return uintx_press(bits, out_bits, in, nin, out, nout); \ -} \ -static int uintx_depress_##bits(uint8_t in_bits, const uint8_t *in, uint64_t nin, \ - uint8_t *out, uint64_t *nout) \ -{ \ - return uintx_depress(in_bits, bits, in, nin, out, nout); \ -} \ - -SLOW5_DEFINE_UINTX(16); -SLOW5_DEFINE_UINTX(32); - -static uint8_t uint_get_minbits(uint64_t max) -{ - uint8_t i; - - for (i = 0; i <= SLOW5_BYTES_TO_BITS(sizeof max); i++) { - if (max < pow(2, i)) - return i; - } - - return SLOW5_BYTES_TO_BITS(sizeof max); -} - -static uint8_t uint_get_minbits_32(const uint32_t *in, uint64_t nin) -{ - uint32_t max; - - max = get_max_u32(in, nin); - return uint_get_minbits(max); -} - -uint64_t uint_bound_32(uint8_t out_bits, uint64_t nin) -{ - uint64_t nin_bytes; - nin_bytes = nin * sizeof (uint32_t); - - return sizeof out_bits + uintx_bound_32(out_bits, nin_bytes); -} - -static int uint_press_32(uint8_t out_bits, const uint32_t *in, uint64_t nin, - uint8_t *out, uint64_t *nout) -{ - int ret; - uint64_t nout_uintx; - const uint8_t *in_uintx; - uint8_t *out_uintx; - - in_uintx = (const uint8_t *) in; - - out[0] = out_bits; - out_uintx = out + 1; - nout_uintx = *nout - sizeof *out; - - ret = uintx_press_32(out_bits, in_uintx, nin, out_uintx, &nout_uintx); - - *nout = 1 + nout_uintx; - return ret; -} - -static int uint_depress_32(const uint8_t *in, uint64_t nin, uint32_t *out, - uint64_t *nout) -{ - int ret; - uint64_t nout_uintx; - const uint8_t *in_uintx; - uint8_t *out_uintx; - uint8_t in_bits; - - in_bits = in[0]; - in_uintx = in + 1; - - out_uintx = (uint8_t *) out; - nout_uintx = *nout * sizeof *out; - - ret = uintx_depress_32(in_bits, in_uintx, nin, out_uintx, &nout_uintx); - - *nout = nout_uintx / sizeof *out; - return ret; -} - -static uint16_t get_max_u16(const uint16_t *in, uint64_t nin) -{ - uint16_t max; - uint64_t i; - - max = 0; - - for (i = 0; i < nin; i++) { - if (in[i] > max) - max = in[i]; - } - - return max; -} - -static uint8_t uint_get_minbits_16(const uint16_t *in, uint64_t nin) -{ - uint16_t max; - - max = get_max_u16(in, nin); - return uint_get_minbits(max); -} - -static uint64_t uint_bound_16(uint8_t out_bits, uint64_t nin) -{ - uint64_t nin_bytes; - nin_bytes = nin * sizeof (uint16_t); - - return sizeof out_bits + uintx_bound_16(out_bits, nin_bytes); -} - -static int uint_press_16(uint8_t out_bits, const uint16_t *in, uint64_t nin, - uint8_t *out, uint64_t *nout) -{ - int ret; - uint64_t nout_uintx; - const uint8_t *in_uintx; - uint8_t *out_uintx; - - in_uintx = (const uint8_t *) in; - - out[0] = out_bits; - out_uintx = out + 1; - nout_uintx = *nout - sizeof *out; - - ret = uintx_press_16(out_bits, in_uintx, nin, out_uintx, &nout_uintx); - - *nout = 1 + nout_uintx; - return ret; -} - -static int uint_depress_16(const uint8_t *in, uint64_t nin, uint16_t *out, - uint64_t *nout) -{ - int ret; - uint64_t nout_uintx; - const uint8_t *in_uintx; - uint8_t *out_uintx; - uint8_t in_bits; - - in_bits = in[0]; - in_uintx = in + 1; - - out_uintx = (uint8_t *) out; - nout_uintx = *nout * sizeof *out; - - ret = uintx_depress_16(in_bits, in_uintx, nin, out_uintx, &nout_uintx); - - *nout = nout_uintx / sizeof *out; - return ret; -} - static void ex_press(const uint16_t *in, uint32_t nin, uint8_t *out, uint64_t *nout) { uint32_t nex; - uint16_t *ex; + uint32_t *ex; uint8_t *ex_press; uint32_t *ex_pos; uint32_t *ex_pos_delta; uint8_t *ex_pos_press; - uint8_t minbits; + //uint8_t minbits; uint64_t nr_press_tmp; uint32_t nex_pos_press; uint32_t nex_press; @@ -1795,10 +1329,13 @@ static void ex_press(const uint16_t *in, uint32_t nin, uint8_t *out, free(ex_pos_press); offset += nex_pos_press; - minbits = uint_get_minbits_16(ex, nex); - nr_press_tmp = uint_bound_16(minbits, nex); + // minbits = uint_get_minbits_16(ex, nex); + // nr_press_tmp = uint_bound_16(minbits, nex); + nr_press_tmp = __slow5_streamvbyte_max_compressedbytes(nex); ex_press = malloc(nr_press_tmp); - (void) uint_press_16(minbits, ex, nex, ex_press, &nr_press_tmp); + //(void) uint_press_16(minbits, ex, nex, ex_press, &nr_press_tmp); + nr_press_tmp=__slow5_streamvbyte_encode(ex, nex, ex_press); + nex_press = (uint32_t) nr_press_tmp; (void) memcpy(out + offset, &nex_press, sizeof nex_press); @@ -1852,14 +1389,14 @@ static void undelta_inplace_increasing_u32(uint32_t *in, uint64_t nin) void ex_depress(const uint8_t *in, uint64_t nin, uint16_t *out, uint32_t *nout) { uint32_t nex; - uint16_t *ex; + uint32_t *ex; uint32_t *ex_pos; uint8_t *ex_pos_press; uint8_t *ex_press; uint32_t nex_press; uint32_t nex_pos_press; - uint64_t nex_pos_delta; - uint64_t nex_cp; + //uint64_t nex_pos_delta; + //uint64_t nex_cp; uint32_t i; uint32_t j; uint64_t offset; @@ -1884,7 +1421,7 @@ void ex_depress(const uint8_t *in, uint64_t nin, uint16_t *out, uint32_t *nout) offset += nex_pos_press; /* TODO don't ignore return */ - nex_pos_delta = nex; + //nex_pos_delta = nex; //fprintf(stderr, "nex_pos_delta = %ld\n", nex_pos_delta); // (void) uint_depress_32(ex_pos_press, nex, ex_pos, &nex_pos_delta); SLOW5_ASSERT(__slow5_streamvbyte_decode(ex_pos_press, ex_pos, nex)==nex_pos_press); @@ -1901,8 +1438,9 @@ void ex_depress(const uint8_t *in, uint64_t nin, uint16_t *out, uint32_t *nout) (void) memcpy(ex_press, in + offset, nex_press); offset += nex_press; - nex_cp = nex; - (void) uint_depress_16(ex_press, nex, ex, &nex_cp); + //nex_cp = nex; + //(void) uint_depress_16(ex_press, nex, ex, &nex_cp); + SLOW5_ASSERT(__slow5_streamvbyte_decode(ex_press, ex, nex)==nex_press); free(ex_press); } else if (nex == 1) { (void) memcpy(ex_pos, in + offset, nex * sizeof *ex_pos);