Merge pull request #92 from sashajenner/ex-zd

qts round

src/slow5_extra.h

@@ -2,6 +2,7 @@
 #define SLOW5_EXTRA_H
 
 #include <dirent.h>
+#include <stdint.h>
 #include <slow5/slow5.h>
 
 #ifdef __cplusplus
@@ -142,6 +143,20 @@ const char *slow5_format_get_str(enum slow5_format format);
 //const uint8_t *str_get_slow5_version(const char *str);
 */
 
+/* Irreversible signal degradation (lossy compression) API */
+
+/*
+ * For array a with length n, zero b LSB of each integer by rounding them to the
+ * nearest power of 2. Set slow5_errno on error.
+ */
+void slow5_arr_qts_round(int16_t *a, uint64_t n, uint8_t b);
+
+/*
+ * For the signal array in a slow5 record, zero b LSB of each integer by
+ * rounding them to the nearest power of 2. Set slow5_errno on error.
+ */
+void slow5_rec_qts_round(struct slow5_rec *r, uint8_t b);
+
 #ifdef __cplusplus
 }
 #endif

src/slow5_press.c

@@ -13,6 +13,7 @@
 #ifdef SLOW5_USE_ZSTD
 #include <zstd.h>
 #endif /* SLOW5_USE_ZSTD */
+#include "slow5_extra.h"
 #include "slow5_misc.h"
 
 extern enum slow5_log_level_opt  slow5_log_level;
@@ -46,6 +47,7 @@ static int16_t *ptr_depress_ex_zd(const uint8_t *ptr, size_t count, size_t *n);
 
 /* other */
 static int vfprintf_compress(struct __slow5_press *comp, FILE *fp, const char *format, va_list ap);
+static int round_to_power_of_2(int number, int number_of_bits);
 
 
 /* Convert the press method back an forth from the library format to spec format
@@ -1956,3 +1958,62 @@ unsigned char *z_inflate_buf(const char *comp_str, size_t *n) {
     return written;
 }
 */
+
+/* Irreversible signal degradation (lossy compression) */
+
+/*
+ * Zero number_of_bits LSB of number by rounding it to the nearest power of 2.
+ * number_of_bits must be >= 1. Return the rounded number.
+ * Taken from https://github.com/hasindu2008/sigtk src/qts.c.
+ */
+static int round_to_power_of_2(int number, int number_of_bits) {
+    //create a binary mask with the specified number of bits
+    int bit_mask = (1 << number_of_bits) - 1;
+
+    //extract out the value of the LSBs considered
+    int lsb_bits = number & bit_mask;
+
+
+    int round_threshold = (1 << (number_of_bits - 1));
+
+    //check if the least significant bits are closer to 0 or 2^n
+    if (lsb_bits < round_threshold) {
+        return (number & ~bit_mask) + 0; //round down to the nearest power of 2
+    } else {
+        return (number & ~bit_mask) + (1 << number_of_bits); //round up to the nearest power of 2
+    }
+}
+
+/*
+ * For array a with length n, zero b LSB of each integer by rounding them to the
+ * nearest power of 2. Set slow5_errno on error.
+ */
+void slow5_arr_qts_round(int16_t *a, uint64_t n, uint8_t b)
+{
+    uint64_t i;
+
+    if (!b)
+        return;
+    if (!a) {
+        SLOW5_ERROR("Argument '%s' cannot be NULL.", SLOW5_TO_STR(a));
+        slow5_errno = SLOW5_ERR_ARG;
+        return;
+    }
+
+    for (i = 0; i < n; i++)
+        a[i] = round_to_power_of_2(a[i], (int) b);
+}
+
+/*
+ * For the signal array in a slow5 record, zero b LSB of each integer by
+ * rounding them to the nearest power of 2. Set slow5_errno on error.
+ */
+void slow5_rec_qts_round(struct slow5_rec *r, uint8_t b)
+{
+    if (!r) {
+        SLOW5_ERROR("Argument '%s' cannot be NULL.", SLOW5_TO_STR(r));
+        slow5_errno = SLOW5_ERR_ARG;
+    } else {
+        slow5_arr_qts_round(r->raw_signal, r->len_raw_signal, b);
+    }
+}