Improvements to gbam C version

.gitignore

@@ -5,3 +5,4 @@ Cargo.lock
 /gbam_tools/gbam_tools/*.so
 /gbam_tools/target
 tests/__pycache__
+*.o

c_attempt/codec_map.json

@@ -0,0 +1,21 @@
+{
+    "COLUMNTYPE_refID": "brotli",
+    "COLUMNTYPE_pos": "brotli",
+    "COLUMNTYPE_mapq": "brotli",
+    "COLUMNTYPE_bin": "brotli",
+    "COLUMNTYPE_flag": "brotli",
+    "COLUMNTYPE_next_refID": "brotli",
+    "COLUMNTYPE_next_pos": "brotli",
+    "COLUMNTYPE_tlen": "brotli",
+    "COLUMNTYPE_read_name": "brotli",
+    "COLUMNTYPE_cigar": "brotli",
+    "COLUMNTYPE_seq": "brotli",
+    "COLUMNTYPE_qual": "brotli",
+    "COLUMNTYPE_tags": "brotli",
+    "COLUMNTYPE_index_read_name": "brotli",
+    "COLUMNTYPE_index_cigar": "brotli",
+    "COLUMNTYPE_index_seq": "brotli",
+    "COLUMNTYPE_index_tags": "brotli",
+    "COLUMNTYPE_index_qual": "brotli"
+  }
+

c_attempt/defs.h

@@ -3,7 +3,7 @@
 #include <stdint.h>
 #include <htslib/sam.h>  
 
-static const int64_t MAX_COLUMN_CHUNK_SIZE = 1024 * 1024 * 150; // 10 MB
+static const int64_t MAX_COLUMN_CHUNK_SIZE = 1024 * 1024 * 6; // 10 MB
 
 static const char* GBAM_MAGIC = "geeBAM20";
 
@@ -58,14 +58,15 @@ typedef struct ColumnChunkMeta
     int64_t file_offset;       // offset in the file where this column starts
     int64_t uncompressed_size; // size of the uncompressed data
     int64_t compressed_size;   // size of the compressed data
-    char codec[5]; // codec used for compression, 0 for no compression
+    char codec[7]; // codec used for compression, 0 for no compression
     struct ColumnChunkMeta *next;      // pointer to the next metadata chunk
     struct ColumnChunkMeta *prev;      // pointer to the next metadata chunk
 } ColumnChunkMeta;
 
 typedef struct Column
 {
     uint8_t *data;
+    int capacity;
     int64_t cur_ptr;
     struct Column* index_column;
 } Column;
@@ -77,6 +78,7 @@ typedef struct  {
     Column *columns;
     ColumnChunkMeta **metadatas; // Pointer to metadata for each column
     int64_t cur_chunk_rec_num[COLUMNTYPE_SIZE] ; // Number of records in the current chunk for each column
+    char codec_map[COLUMNTYPE_SIZE][8];
 } Writer;
 
 typedef struct {
@@ -101,6 +103,7 @@ struct bam_hdr_t ;
 
 Writer *create_writer(FILE* fd, bam_hdr_t  *header);
 int write_bam_record(Writer *writer, bam1_t *aln);
+void ensure_column_capacity(Column* col, int needed_size);
 void close_writer(Writer *writer);
 
 Reader *make_reader(char *fp);

c_attempt/main.c

@@ -7,13 +7,18 @@
 #include <unistd.h>
 #include <stdio.h>
 #include <sys/stat.h>
+#include <time.h>
+
+void print_progress_eta(int count, time_t start_time);
 
 int write_gbam(char* file_path){
+    time_t start_time = time(NULL);
+    int record_count = 0;
     // Open the file for writing
     FILE *fp = fopen(file_path, "wb");
     if (!fp) {
         perror("Failed to open file for writing");
-        return;
+        return 1;
     }
 
     samFile *in = sam_open("-", "r");
@@ -24,7 +29,7 @@ int write_gbam(char* file_path){
 
     if (hts_set_threads(in, 8) != 0) {
         sam_close(fp);
-        return NULL;
+        return -1;
     }
 
     bam_hdr_t *header = sam_hdr_read(in);
@@ -62,8 +67,12 @@ int write_gbam(char* file_path){
             fclose(fp);
             return 1;
         }
+        record_count++;
+        print_progress_eta(record_count, start_time);
     }
 
+    printf("\nCompleted writing %d records.\n", record_count);
+
     close_writer(writer);
 
     // Cleanup
@@ -106,6 +115,10 @@ int read_gbam(char* file_path){
     setvbuf(stdout, stdout_buffer, _IOFBF, sizeof(stdout_buffer));
     kstring_t str = {0, 0, NULL};
 
+    htsFile *fp = hts_open("-", "w");  // "-" means stdout
+    sam_hdr_write(fp, reader->header);
+    hts_close(fp);
+
     for (int i = 0; i < reader->rec_num; i++) {
         read_record(reader, i, aln);
 
@@ -126,6 +139,25 @@ int read_gbam(char* file_path){
     return 0;
 }
 
+void print_progress_eta(int count, time_t start_time) {
+    if (count % 1000 != 0) return;  // throttle updates
+
+    time_t now = time(NULL);
+    double elapsed = difftime(now, start_time);
+    if (elapsed < 1.0) elapsed = 1.0;  // avoid division by zero
+
+    double records_per_sec = count / elapsed;
+    int est_total = (int)(elapsed > 5 ? (count * 1.3) : (count * 2));  // dynamic guess
+    int est_remaining = est_total - count;
+    int eta_seconds = (int)(est_remaining / records_per_sec);
+
+    int mins = eta_seconds / 60;
+    int secs = eta_seconds % 60;
+
+    printf("\rProcessed: %d | Speed: %.1f rec/s | ETA: %02d:%02d", count, records_per_sec, mins, secs);
+    fflush(stdout);
+}
+
 int main(int argc, char *argv[]) {
 
     if (argc < 2) {

c_attempt/makefile

@@ -1,8 +1,11 @@
 # Compiler and flags
 CC = gcc
-CFLAGS = -O3 -g -fno-omit-frame-pointer 
+# CFLAGS = -O3 -g -fno-omit-frame-pointer -fsanitize=address
 
-LDFLAGS = -ljson-c -lhts -lz -llz4
+# LDFLAGS = -ljson-c -lhts -lz -llz4 -lbrotlidec -lbrotlienc -lzstd
+
+CFLAGS = -O3 -g -fno-omit-frame-pointer -I$(HOME)/.local/include -fopenmp
+LDFLAGS = -L$(HOME)/.local/lib -ljson-c -lhts -lz -llz4 -lbrotlidec -lbrotlienc -lbrotlicommon -lzstd -fopenmp
 
 # Source and object files
 SRCS = main.c reader.c writer.c

c_attempt/meta.h

@@ -45,10 +45,11 @@ void write_meta(FILE *fp, ColumnChunkMeta **array, size_t size) {
     fprintf(fp, "}\n\0");
 }
 
-void write_header(FILE *fp, int64_t seekpos, int64_t meta_size) {
+void write_header(FILE *fp, int64_t seekpos, int64_t meta_size, int64_t header_len) {
     fprintf(fp, "{");
     fprintf(fp, "\"seekpos\": %lld,\n", (long long)seekpos);
-    fprintf(fp, "\"meta_size\": %lld\n", (long long)meta_size);
+    fprintf(fp, "\"meta_size\": %lld,\n", (long long)meta_size);
+    fprintf(fp, "\"header_len\": %lld\n", (long long)header_len);
     fprintf(fp, "}\0");
 }
 
@@ -99,5 +100,5 @@ void parse_meta_from_json_string(char *json_str, Reader *reader) {
 
     reader->metadatas = array;
 
-    return array;
+    // return array;
 }

c_attempt/reader.c

@@ -5,13 +5,18 @@
 #include <htslib/hts.h>
 #include <stdbool.h>
 #include <zlib.h>
+#include <lz4.h>
+#include <brotli/decode.h>
+#include <zstd.h>
 
 
 #define ADJUSTED_OFFSET(COLUMNTYPE) \
     (rec_num-(reader->loaded_since_rec_num[COLUMNTYPE]))
 
 #define bam_reg2bin(beg,end) hts_reg2bin((beg),(end),14,5)
 
+void parse_meta_from_json_string(const char *json_str, Reader *reader);
+
 
 static void bam_cigar2rqlens(int n_cigar, const uint32_t *cigar,
                              hts_pos_t *rlen, hts_pos_t *qlen)
@@ -47,6 +52,13 @@ Reader* make_reader(char* file){
         meta_size = json_object_get_int64(meta_size_obj);
     }
 
+    int64_t header_len = -1;
+    struct json_object *header_len_obj;
+    if (json_object_object_get_ex(root, "header_len", &header_len_obj)) {
+        header_len = json_object_get_int64(header_len_obj);
+    }
+    assert(header_len > 0);
+
     json_object_put(root); 
 
 
@@ -58,7 +70,12 @@ Reader* make_reader(char* file){
     parse_meta_from_json_string(file+seekpos, reader);
 
     reader->mmaped_file = file;
-    reader->header = sam_hdr_parse(strlen(file+seekpos+meta_size), file+seekpos+meta_size);
+    // reader->header = sam_hdr_parse(strlen(file+seekpos+meta_size), file+seekpos+meta_size); <- this doesn't work
+    // because strlen(...) calculates length until the first \0 byte, which is not guaranteed in binary data.
+    int32_t *header_len_ptr = (int32_t *)(file + seekpos + meta_size);
+    char *header_start = (char *)(header_len_ptr + 1);
+    assert(*header_len_ptr == header_len);  // Sanity check to confirm the calculated header length is matched with the header length stored in the metadata.
+    reader->header = sam_hdr_parse(*header_len_ptr, header_start);
     reader->rec_num = 0;
     reader->columns = (Column*)calloc(COLUMNTYPE_SIZE, sizeof(Column));
 
@@ -147,7 +164,36 @@ void fetch_field(Reader* reader, int64_t rec_num, int64_t COLUMNTYPE){
                                     meta->uncompressed_size);
             if (decompressed_size < 0) {
                 printf("LZ4 decompression failed with error code: %d\n", decompressed_size);
-                return -1;
+                return;
+            }
+        }
+        else if(strcmp(meta->codec, "brotli") == 0){
+            size_t decompressed_size = meta->uncompressed_size;
+            size_t compressed_size = meta->compressed_size;
+
+            if (reader->m_chunk_memory[COLUMNTYPE] < decompressed_size) {
+                // resize buffer to expected size
+                free(reader->columns[COLUMNTYPE].data);
+                reader->columns[COLUMNTYPE].data = (uint8_t*)malloc(decompressed_size);
+                reader->m_chunk_memory[COLUMNTYPE] = decompressed_size;
+            }
+
+            BrotliDecoderResult res = BrotliDecoderDecompress(
+                compressed_size, (const uint8_t*)read_buffer,
+                &decompressed_size, reader->columns[COLUMNTYPE].data);
+
+            if (res != BROTLI_DECODER_RESULT_SUCCESS || decompressed_size != meta->uncompressed_size) {
+                fprintf(stderr, "Failed to decompress Brotli data or unexpected size (got %zu, expected %zu)\n",
+                        decompressed_size, meta->uncompressed_size);
+                exit(1);
+            }
+        }
+        else if(strcmp(meta->codec, "zstd") == 0){
+            size_t result = ZSTD_decompress(reader->columns[COLUMNTYPE].data, meta->uncompressed_size,
+                                            read_buffer, meta->compressed_size);
+            if (ZSTD_isError(result) || result != meta->uncompressed_size) {
+                fprintf(stderr, "Zstd decompression failed: %s\n", ZSTD_getErrorName(result));
+                exit(1);
             }
         }
         else{
@@ -243,10 +289,10 @@ void read_record(Reader* reader, int64_t rec_num, bam1_t* aln) {
     {
         default_name = true;
         l_qname = 1;
-        qname_nuls = 4 - l_qname % 4;
+        qname_nuls = 0;
     }
     else{
-        qname_nuls = 4 - l_qname % 4;
+        qname_nuls = 0;
     }
 
     uint64_t bytes_we_need = l_qname + qname_nuls + l_cigar + l_seq + l_qual + l_tags;
@@ -268,9 +314,7 @@ void read_record(Reader* reader, int64_t rec_num, bam1_t* aln) {
     }
     else{
         memcpy(aln->data, 
-           &reader->columns[COLUMNTYPE_read_name].data[read_name_beg], l_qname);
-        for (int i = 0; i < qname_nuls; i++) 
-            aln->data[l_qname + i] = '\0'; // Fill with null bytes
+            &reader->columns[COLUMNTYPE_read_name].data[read_name_beg], l_qname);
     }
 
     memcpy(&aln->data[l_qname + qname_nuls], 
@@ -294,16 +338,14 @@ void read_record(Reader* reader, int64_t rec_num, bam1_t* aln) {
     aln->core.tid = refID;
     aln->core.bin = bam_reg2bin(pos, pos + rlen);
     aln->core.qual = mapq;
-    aln->core.l_extranul = (uint8_t)(qname_nuls - 1);
+    aln->core.l_extranul = 0;
     aln->core.flag = flag;
-    aln->core.l_qname = (uint16_t)(l_qname + qname_nuls);
+    aln->core.l_qname = (uint16_t)(l_qname);
     aln->core.n_cigar = (uint32_t)(l_cigar >> 2); // l_cigar is in bytes, n_cigar is in 32-bit words
     aln->core.l_qseq = (int32_t)l_qual;
     aln->core.mtid = next_refID;
     aln->core.mpos = next_pos;
     aln->core.isize = tlen;
-
-    return aln;
 }
 
 void close_reader(Reader *reader) {