VCF-42: Reduce memory consumption when a single position returns many samples (#867)

* Added external MD5 library

* Added a join utility function

This function joins a vector of strings using a given delimiter character while optionally omitting empty strings.

* Added a SampleHeaders class to abstract the representation of loaded sample headers

SampleHeaders is a subclass of the TileDBVCFDataset class since it is intended to be instantiated and loaded exclusively by TileDBVCFDataset. Internally SampleHeaders uses MD5 checksums and maps to only load unique headers and to associate samples with these unique headers, respectively.

* Updated relevant code paths to use new SampleHeaders class

* Updated VCF store unit tests to use new SampleHeaders class

* unique_ptr instances returned by the SampleHeaders class are no longer released

This is to prevent memory leaks while making the lifetime of these managed pointers more obvious.

* Added a samples_view() method to the SampleHeaders class

This allows the names of samples stored in a SampleHeaders instance to be efficiently iterated via a view of the keys of an internal map.

* Updated Reader to use samples_view() method of SampleHeaders instance

This allows the names of only samples that have headers loaded to be iterated efficiently, both in terms of run-time and memory consumption.

* Fixed MD5 stack overflow bug on macOS

Specifically, the std::string(const char* s, size_t n) constructor on macOS seems to ignore the n parameter and overflows. This was mitigated by adding a terminating character to s and using the std::string(const char* s) constructor instead.

Author: alancleary

libtiledbvcf/external/md5/md5.cc

@@ -0,0 +1,215 @@
+#include "md5.h"
+#include <cstddef>
+
+namespace md5 {
+
+/*
+ * Constants defined by the MD5 algorithm
+ */
+#define A 0x67452301
+#define B 0xefcdab89
+#define C 0x98badcfe
+#define D 0x10325476
+
+static uint32_t S[] = {7,  12, 17, 22, 7,  12, 17, 22, 7,  12, 17, 22, 7,
+                       12, 17, 22, 5,  9,  14, 20, 5,  9,  14, 20, 5,  9,
+                       14, 20, 5,  9,  14, 20, 4,  11, 16, 23, 4,  11, 16,
+                       23, 4,  11, 16, 23, 4,  11, 16, 23, 6,  10, 15, 21,
+                       6,  10, 15, 21, 6,  10, 15, 21, 6,  10, 15, 21};
+
+static uint32_t K[] = {
+    0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf, 0x4787c62a,
+    0xa8304613, 0xfd469501, 0x698098d8, 0x8b44f7af, 0xffff5bb1, 0x895cd7be,
+    0x6b901122, 0xfd987193, 0xa679438e, 0x49b40821, 0xf61e2562, 0xc040b340,
+    0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x02441453, 0xd8a1e681, 0xe7d3fbc8,
+    0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905, 0xfcefa3f8,
+    0x676f02d9, 0x8d2a4c8a, 0xfffa3942, 0x8771f681, 0x6d9d6122, 0xfde5380c,
+    0xa4beea44, 0x4bdecfa9, 0xf6bb4b60, 0xbebfbc70, 0x289b7ec6, 0xeaa127fa,
+    0xd4ef3085, 0x04881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8, 0xc4ac5665,
+    0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3, 0x8f0ccc92,
+    0xffeff47d, 0x85845dd1, 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314, 0x4e0811a1,
+    0xf7537e82, 0xbd3af235, 0x2ad7d2bb, 0xeb86d391};
+
+/*
+ * Padding used to make the size (in bits) of the input congruent to 448 mod 512
+ */
+static uint8_t PADDING[] = {
+    0x80, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+    0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00};
+
+/*
+ * Bit-manipulation functions defined by the MD5 algorithm
+ */
+#define F(X, Y, Z) ((X & Y) | (~X & Z))
+#define G(X, Y, Z) ((X & Z) | (Y & ~Z))
+#define H(X, Y, Z) (X ^ Y ^ Z)
+#define I(X, Y, Z) (Y ^ (X | ~Z))
+
+/*
+ * Rotates a 32-bit word left by n bits
+ */
+uint32_t rotateLeft(uint32_t x, uint32_t n) {
+  return (x << n) | (x >> (32 - n));
+}
+
+void md5_init(MD5Context* ctx) {
+  ctx->size = (uint64_t)0;
+
+  ctx->buffer[0] = (uint32_t)A;
+  ctx->buffer[1] = (uint32_t)B;
+  ctx->buffer[2] = (uint32_t)C;
+  ctx->buffer[3] = (uint32_t)D;
+}
+
+void md5_update(
+    MD5Context* ctx, const uint8_t* input_buffer, size_t input_len) {
+  uint32_t input[16];
+  unsigned int offset = ctx->size % 64;
+  ctx->size += (uint64_t)input_len;
+
+  // Copy each byte in input_buffer into the next space in our context input
+  for (unsigned int i = 0; i < input_len; ++i) {
+    ctx->input[offset++] = (uint8_t)*(input_buffer + i);
+
+    // If we've filled our context input, copy it into our local array input
+    // then reset the offset to 0 and fill in a new buffer.
+    // Every time we fill out a chunk, we run it through the algorithm
+    // to enable some back and forth between cpu and i/o
+    if (offset % 64 == 0) {
+      for (unsigned int j = 0; j < 16; ++j) {
+        // Convert to little-endian
+        // The local variable `input` our 512-bit chunk separated into 32-bit
+        // words we can use in calculations
+        input[j] = (uint32_t)(ctx->input[(j * 4) + 3]) << 24 |
+                   (uint32_t)(ctx->input[(j * 4) + 2]) << 16 |
+                   (uint32_t)(ctx->input[(j * 4) + 1]) << 8 |
+                   (uint32_t)(ctx->input[(j * 4)]);
+      }
+      md5_step(ctx->buffer, input);
+      offset = 0;
+    }
+  }
+}
+
+void md5_finalize(MD5Context* ctx) {
+  uint32_t input[16];
+  unsigned int offset = ctx->size % 64;
+  unsigned int padding_length = offset < 56 ? 56 - offset : (56 + 64) - offset;
+
+  // Fill in the padding and undo the changes to size that resulted from the
+  // update
+  md5_update(ctx, PADDING, padding_length);
+  ctx->size -= (uint64_t)padding_length;
+
+  // Do a final update (internal to this function)
+  // Last two 32-bit words are the two halves of the size (converted from bytes
+  // to bits)
+  for (unsigned int j = 0; j < 14; ++j) {
+    input[j] = (uint32_t)(ctx->input[(j * 4) + 3]) << 24 |
+               (uint32_t)(ctx->input[(j * 4) + 2]) << 16 |
+               (uint32_t)(ctx->input[(j * 4) + 1]) << 8 |
+               (uint32_t)(ctx->input[(j * 4)]);
+  }
+  input[14] = (uint32_t)(ctx->size * 8);
+  input[15] = (uint32_t)((ctx->size * 8) >> 32);
+
+  md5_step(ctx->buffer, input);
+
+  // Move the result into digest (convert from little-endian)
+  for (unsigned int i = 0; i < 4; ++i) {
+    ctx->digest[(i * 4) + 0] = (uint8_t)((ctx->buffer[i] & 0x000000FF));
+    ctx->digest[(i * 4) + 1] = (uint8_t)((ctx->buffer[i] & 0x0000FF00) >> 8);
+    ctx->digest[(i * 4) + 2] = (uint8_t)((ctx->buffer[i] & 0x00FF0000) >> 16);
+    ctx->digest[(i * 4) + 3] = (uint8_t)((ctx->buffer[i] & 0xFF000000) >> 24);
+  }
+}
+
+void md5_step(uint32_t* buffer, uint32_t* input) {
+  uint32_t AA = buffer[0];
+  uint32_t BB = buffer[1];
+  uint32_t CC = buffer[2];
+  uint32_t DD = buffer[3];
+
+  uint32_t E;
+
+  unsigned int j;
+
+  for (unsigned int i = 0; i < 64; ++i) {
+    switch (i / 16) {
+      case 0:
+        E = F(BB, CC, DD);
+        j = i;
+        break;
+      case 1:
+        E = G(BB, CC, DD);
+        j = ((i * 5) + 1) % 16;
+        break;
+      case 2:
+        E = H(BB, CC, DD);
+        j = ((i * 3) + 5) % 16;
+        break;
+      default:
+        E = I(BB, CC, DD);
+        j = (i * 7) % 16;
+        break;
+    }
+
+    uint32_t temp = DD;
+    DD = CC;
+    CC = BB;
+    BB = BB + rotateLeft(AA + E + K[i] + input[j], S[i]);
+    AA = temp;
+  }
+
+  buffer[0] += AA;
+  buffer[1] += BB;
+  buffer[2] += CC;
+  buffer[3] += DD;
+}
+
+void md5_string(const char* input, uint8_t* result) {
+  MD5Context ctx;
+  md5_init(&ctx);
+  md5_update(&ctx, (uint8_t*)input, strlen(input));
+  md5_finalize(&ctx);
+
+  memcpy(result, ctx.digest, 16);
+}
+
+void md5_file(FILE* file, uint8_t* result) {
+  char* input_buffer = (char*)malloc(1024);
+  size_t input_size = 0;
+
+  MD5Context ctx;
+  md5_init(&ctx);
+
+  while ((input_size = fread(input_buffer, 1, 1024, file)) > 0) {
+    md5_update(&ctx, (uint8_t*)input_buffer, input_size);
+  }
+
+  md5_finalize(&ctx);
+
+  free(input_buffer);
+
+  memcpy(result, ctx.digest, 16);
+}
+
+void md5_to_hex(const uint8_t* input, char* output) {
+  for (size_t i = 0; i < 16; i++) {
+    sprintf(output, "%02x", input[i]);
+    output += 2;
+  }
+}
+
+std::string md5_to_hex(const uint8_t* input) {
+  char output[33];
+  md5_to_hex(input, output);
+  output[32] = '\0';
+  return std::string(output);
+}
+
+}  // namespace md5

libtiledbvcf/external/md5/md5.h

@@ -0,0 +1,77 @@
+/*
+ * This file (md5.h) and the corresponding md5.cc file are licensed under the
+ * Unlicense license and are therefore released into the public domain.
+ *
+ * This is an a adaptation of the MD5 reference implementation from
+ * https://github.com/Zunawe/md5-c. It is derived from the RSA Data Security,
+ * Inc. MD5 Message-Digest Algorithm and modified slightly to be functionally
+ * identical but condensed into control structures.
+ */
+
+#ifndef TILEDB_EXTERNAL_MD5_H
+#define TILEDB_EXTERNAL_MD5_H
+
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <string>
+
+namespace md5 {
+
+typedef struct {
+  uint64_t size;       // Size of input in bytes
+  uint32_t buffer[4];  // Current accumulation of hash
+  uint8_t input[64];   // Input to be used in the next step
+  uint8_t digest[16];  // Result of algorithm
+} MD5Context;
+
+/*
+ * Initialize an md5 context
+ */
+void md5_init(MD5Context* ctx);
+
+/*
+ * Add some amount of input to the context
+ *
+ * If the input fills out a block of 512 bits, apply the algorithm (md5_step)
+ * and save the result in the buffer. Also updates the overall size.
+ */
+void md5_update(MD5Context* ctx, const uint8_t* input, size_t input_len);
+
+/*
+ * Pad the current input to get to 448 bytes, append the size in bits to the
+ * very end, and save the result of the final iteration into digest.
+ */
+void md5_finalize(MD5Context* ctx);
+
+/*
+ * Step on 512 bits of input with the main MD5 algorithm.
+ */
+void md5_step(uint32_t* buffer, uint32_t* input);
+
+/*
+ * Functions that run the algorithm on the provided input and put the digest
+ * into result. result should be able to store 16 bytes.
+ */
+void md5_string(const char* input, uint8_t* result);
+
+/*
+ * Functions that run the algorithm on the contents of the provided file and
+ * put the digest into result. result should be able to store 16 bytes.
+ */
+void md5_file(FILE* file, uint8_t* result);
+
+/*
+ * Converts an md5 array into a C string.
+ */
+void md5_to_hex(const uint8_t* input, char* output);
+
+/*
+ * Converts an md5 array into a string.
+ */
+std::string md5_to_hex(const uint8_t* input);
+
+}  // namespace md5
+
+#endif

libtiledbvcf/src/CMakeLists.txt

@@ -90,6 +90,7 @@ set(TILEDB_VCF_SOURCES
 
 set(TILEDB_VCF_EXTERNAL_SOURCES
   ${CMAKE_CURRENT_SOURCE_DIR}/../external/base64/base64.cc
+  ${CMAKE_CURRENT_SOURCE_DIR}/../external/md5/md5.cc
 )
 
 add_library(TILEDB_VCF_OBJECTS OBJECT