bgzf.c

   1 /*
   2  * The Broad Institute
   3  * SOFTWARE COPYRIGHT NOTICE AGREEMENT
   4  * This software and its documentation are copyright 2008 by the
   5  * Broad Institute/Massachusetts Institute of Technology. All rights are reserved.
   6  *
   7  * This software is supplied without any warranty or guaranteed support whatsoever.
   8  * Neither the Broad Institute nor MIT can be responsible for its use, misuse,
   9  * or functionality.
  10  */
  11
  12 /* 2009-06-12 by lh3: support a mode string like "wu" where 'u' for uncompressed output */
  13
  14 #include <stdio.h>
  15 #include <stdlib.h>
  16 #include <string.h>
  17 #include <unistd.h>
  18 #include <fcntl.h>
  19 #include <sys/types.h>
  20 #include <sys/stat.h>
  21 #include "bgzf.h"
  22
  23 extern off_t ftello(FILE *stream);
  24 extern int fseeko(FILE *stream, off_t offset, int whence);
  25
  26 typedef int8_t byte;
  27
  28 static const int DEFAULT_BLOCK_SIZE = 64 * 1024;
  29 static const int MAX_BLOCK_SIZE = 64 * 1024;
  30
  31 static const int BLOCK_HEADER_LENGTH = 18;
  32 static const int BLOCK_FOOTER_LENGTH = 8;
  33
  34 static const int GZIP_ID1 = 31;
  35 static const int GZIP_ID2 = 139;
  36 static const int CM_DEFLATE = 8;
  37 static const int FLG_FEXTRA = 4;
  38 static const int OS_UNKNOWN = 255;
  39 static const int BGZF_ID1 = 66; // 'B'
  40 static const int BGZF_ID2 = 67; // 'C'
  41 static const int BGZF_LEN = 2;
  42 static const int BGZF_XLEN = 6; // BGZF_LEN+4
  43
  44 static const int GZIP_WINDOW_BITS = -15; // no zlib header
  45 static const int Z_DEFAULT_MEM_LEVEL = 8;
  46
  47
  48 inline
  49 void
  50 packInt16(uint8_t* buffer, uint16_t value)
  51 {
  52     buffer[0] = value;
  53     buffer[1] = value >> 8;
  54 }
  55
  56 inline
  57 int
  58 unpackInt16(const uint8_t* buffer)
  59 {
  60     return (buffer[0] | (buffer[1] << 8));
  61 }
  62
  63 inline
  64 void
  65 packInt32(uint8_t* buffer, uint32_t value)
  66 {
  67     buffer[0] = value;
  68     buffer[1] = value >> 8;
  69     buffer[2] = value >> 16;
  70     buffer[3] = value >> 24;
  71 }
  72
  73 inline
  74 int
  75 min(int x, int y)
  76 {
  77     return (x < y) ? x : y;
  78 }
  79
  80 static
  81 void
  82 report_error(BGZF* fp, const char* message) {
  83     fp->error = message;
  84 }
  85
  86 static BGZF *bgzf_read_init()
  87 {
  88         BGZF *fp;
  89         fp = calloc(1, sizeof(BGZF));
  90     fp->uncompressed_block_size = MAX_BLOCK_SIZE;
  91     fp->uncompressed_block = malloc(MAX_BLOCK_SIZE);
  92     fp->compressed_block_size = MAX_BLOCK_SIZE;
  93     fp->compressed_block = malloc(MAX_BLOCK_SIZE);
  94         return fp;
  95 }
  96
  97 static
  98 BGZF*
  99 open_read(int fd)
 100 {
 101 #ifdef _USE_KNETFILE
 102     knetFile *file = knet_dopen(fd, "r");
 103 #else
 104     FILE* file = fdopen(fd, "r");
 105 #endif
 106     BGZF* fp;
 107         if (file == 0) return 0;
 108         fp = bgzf_read_init();
 109     fp->file_descriptor = fd;
 110     fp->open_mode = 'r';
 111 #ifdef _USE_KNETFILE
 112     fp->x.fpr = file;
 113 #else
 114     fp->file = file;
 115 #endif
 116     return fp;
 117 }
 118
 119 static
 120 BGZF*
 121 open_write(int fd, bool is_uncompressed)
 122 {
 123     FILE* file = fdopen(fd, "w");
 124     BGZF* fp;
 125         if (file == 0) return 0;
 126         fp = malloc(sizeof(BGZF));
 127     fp->file_descriptor = fd;
 128     fp->open_mode = 'w';
 129     fp->owned_file = 0; fp->is_uncompressed = is_uncompressed;
 130 #ifdef _USE_KNETFILE
 131     fp->x.fpw = file;
 132 #else
 133     fp->file = file;
 134 #endif
 135     fp->uncompressed_block_size = DEFAULT_BLOCK_SIZE;
 136     fp->uncompressed_block = NULL;
 137     fp->compressed_block_size = MAX_BLOCK_SIZE;
 138     fp->compressed_block = malloc(MAX_BLOCK_SIZE);
 139     fp->block_address = 0;
 140     fp->block_offset = 0;
 141     fp->block_length = 0;
 142     fp->error = NULL;
 143     return fp;
 144 }
 145
 146 BGZF*
 147 bgzf_open(const char* __restrict path, const char* __restrict mode)
 148 {
 149     BGZF* fp = NULL;
 150     if (mode[0] == 'r' || mode[0] == 'R') { /* The reading mode is preferred. */
 151 #ifdef _USE_KNETFILE
 152                 knetFile *file = knet_open(path, mode);
 153                 if (file == 0) return 0;
 154                 fp = bgzf_read_init();
 155                 fp->file_descriptor = -1;
 156                 fp->open_mode = 'r';
 157                 fp->x.fpr = file;
 158 #else
 159                 int oflag = O_RDONLY;
 160                 int fd = open(path, oflag);
 161                 if (fd == -1) return 0;
 162         fp = open_read(fd);
 163 #endif
 164     } else if (mode[0] == 'w' || mode[0] == 'W') {
 165                 int oflag = O_WRONLY | O_CREAT | O_TRUNC;
 166                 int fd = open(path, oflag, 0644);
 167                 if (fd == -1) return 0;
 168         fp = open_write(fd, strstr(mode, "u")? 1 : 0);
 169     }
 170     if (fp != NULL) {
 171         fp->owned_file = 1;
 172     }
 173     return fp;
 174 }
 175
 176 BGZF*
 177 bgzf_fdopen(int fd, const char * __restrict mode)
 178 {
 179         if (fd == -1) return 0;
 180     if (mode[0] == 'r' || mode[0] == 'R') {
 181         return open_read(fd);
 182     } else if (mode[0] == 'w' || mode[0] == 'W') {
 183         return open_write(fd, strstr(mode, "u")? 1 : 0);
 184     } else {
 185         return NULL;
 186     }
 187 }
 188
 189 static
 190 int
 191 deflate_block(BGZF* fp, int block_length)
 192 {
 193     // Deflate the block in fp->uncompressed_block into fp->compressed_block.
 194     // Also adds an extra field that stores the compressed block length.
 195
 196     byte* buffer = fp->compressed_block;
 197     int buffer_size = fp->compressed_block_size;
 198
 199     // Init gzip header
 200     buffer[0] = GZIP_ID1;
 201     buffer[1] = GZIP_ID2;
 202     buffer[2] = CM_DEFLATE;
 203     buffer[3] = FLG_FEXTRA;
 204     buffer[4] = 0; // mtime
 205     buffer[5] = 0;
 206     buffer[6] = 0;
 207     buffer[7] = 0;
 208     buffer[8] = 0;
 209     buffer[9] = OS_UNKNOWN;
 210     buffer[10] = BGZF_XLEN;
 211     buffer[11] = 0;
 212     buffer[12] = BGZF_ID1;
 213     buffer[13] = BGZF_ID2;
 214     buffer[14] = BGZF_LEN;
 215     buffer[15] = 0;
 216     buffer[16] = 0; // placeholder for block length
 217     buffer[17] = 0;
 218
 219     // loop to retry for blocks that do not compress enough
 220     int input_length = block_length;
 221     int compressed_length = 0;
 222     while (1) {
 223                 int compress_level = fp->is_uncompressed? 0 : Z_DEFAULT_COMPRESSION;
 224         z_stream zs;
 225         zs.zalloc = NULL;
 226         zs.zfree = NULL;
 227         zs.next_in = fp->uncompressed_block;
 228         zs.avail_in = input_length;
 229         zs.next_out = (void*)&buffer[BLOCK_HEADER_LENGTH];
 230         zs.avail_out = buffer_size - BLOCK_HEADER_LENGTH - BLOCK_FOOTER_LENGTH;
 231
 232         int status = deflateInit2(&zs, compress_level, Z_DEFLATED,
 233                                   GZIP_WINDOW_BITS, Z_DEFAULT_MEM_LEVEL, Z_DEFAULT_STRATEGY);
 234         if (status != Z_OK) {
 235             report_error(fp, "deflate init failed");
 236             return -1;
 237         }
 238         status = deflate(&zs, Z_FINISH);
 239         if (status != Z_STREAM_END) {
 240             deflateEnd(&zs);
 241             if (status == Z_OK) {
 242                 // Not enough space in buffer.
 243                 // Can happen in the rare case the input doesn't compress enough.
 244                 // Reduce the amount of input until it fits.
 245                 input_length -= 1024;
 246                 if (input_length <= 0) {
 247                     // should never happen
 248                     report_error(fp, "input reduction failed");
 249                     return -1;
 250                 }
 251                 continue;
 252             }
 253             report_error(fp, "deflate failed");
 254             return -1;
 255         }
 256         status = deflateEnd(&zs);
 257         if (status != Z_OK) {
 258             report_error(fp, "deflate end failed");
 259             return -1;
 260         }
 261         compressed_length = zs.total_out;
 262         compressed_length += BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH;
 263         if (compressed_length > MAX_BLOCK_SIZE) {
 264             // should never happen
 265             report_error(fp, "deflate overflow");
 266             return -1;
 267         }
 268         break;
 269     }
 270
 271     packInt16((uint8_t*)&buffer[16], compressed_length-1);
 272     uint32_t crc = crc32(0L, NULL, 0L);
 273     crc = crc32(crc, fp->uncompressed_block, input_length);
 274     packInt32((uint8_t*)&buffer[compressed_length-8], crc);
 275     packInt32((uint8_t*)&buffer[compressed_length-4], input_length);
 276
 277     int remaining = block_length - input_length;
 278     if (remaining > 0) {
 279         if (remaining > input_length) {
 280             // should never happen (check so we can use memcpy)
 281             report_error(fp, "remainder too large");
 282             return -1;
 283         }
 284         memcpy(fp->uncompressed_block,
 285                fp->uncompressed_block + input_length,
 286                remaining);
 287     }
 288     fp->block_offset = remaining;
 289     return compressed_length;
 290 }
 291
 292 static
 293 int
 294 inflate_block(BGZF* fp, int block_length)
 295 {
 296     // Inflate the block in fp->compressed_block into fp->uncompressed_block
 297
 298     z_stream zs;
 299     zs.zalloc = NULL;
 300     zs.zfree = NULL;
 301     zs.next_in = fp->compressed_block + 18;
 302     zs.avail_in = block_length - 16;
 303     zs.next_out = fp->uncompressed_block;
 304     zs.avail_out = fp->uncompressed_block_size;
 305
 306     int status = inflateInit2(&zs, GZIP_WINDOW_BITS);
 307     if (status != Z_OK) {
 308         report_error(fp, "inflate init failed");
 309         return -1;
 310     }
 311     status = inflate(&zs, Z_FINISH);
 312     if (status != Z_STREAM_END) {
 313         inflateEnd(&zs);
 314         report_error(fp, "inflate failed");
 315         return -1;
 316     }
 317     status = inflateEnd(&zs);
 318     if (status != Z_OK) {
 319         report_error(fp, "inflate failed");
 320         return -1;
 321     }
 322     return zs.total_out;
 323 }
 324
 325 static
 326 int
 327 check_header(const byte* header)
 328 {
 329     return (header[0] == GZIP_ID1 &&
 330             header[1] == (byte) GZIP_ID2 &&
 331             header[2] == Z_DEFLATED &&
 332             (header[3] & FLG_FEXTRA) != 0 &&
 333             unpackInt16((uint8_t*)&header[10]) == BGZF_XLEN &&
 334             header[12] == BGZF_ID1 &&
 335             header[13] == BGZF_ID2 &&
 336             unpackInt16((uint8_t*)&header[14]) == BGZF_LEN);
 337 }
 338
 339 static
 340 int
 341 read_block(BGZF* fp)
 342 {
 343     byte header[BLOCK_HEADER_LENGTH];
 344 #ifdef _USE_KNETFILE
 345     int64_t block_address = knet_tell(fp->x.fpr);
 346     int count = knet_read(fp->x.fpr, header, sizeof(header));
 347 #else
 348     int64_t block_address = ftello(fp->file);
 349     int count = fread(header, 1, sizeof(header), fp->file);
 350 #endif
 351     if (count == 0) {
 352         fp->block_length = 0;
 353         return 0;
 354     }
 355     if (count != sizeof(header)) {
 356         report_error(fp, "read failed");
 357         return -1;
 358     }
 359     if (!check_header(header)) {
 360         report_error(fp, "invalid block header");
 361         return -1;
 362     }
 363     int block_length = unpackInt16((uint8_t*)&header[16]) + 1;
 364     byte* compressed_block = (byte*) fp->compressed_block;
 365     memcpy(compressed_block, header, BLOCK_HEADER_LENGTH);
 366     int remaining = block_length - BLOCK_HEADER_LENGTH;
 367 #ifdef _USE_KNETFILE
 368     count = knet_read(fp->x.fpr, &compressed_block[BLOCK_HEADER_LENGTH], remaining);
 369 #else
 370     count = fread(&compressed_block[BLOCK_HEADER_LENGTH], 1, remaining, fp->file);
 371 #endif
 372     if (count != remaining) {
 373         report_error(fp, "read failed");
 374         return -1;
 375     }
 376     count = inflate_block(fp, block_length);
 377     if (count < 0) {
 378         return -1;
 379     }
 380     if (fp->block_length != 0) {
 381         // Do not reset offset if this read follows a seek.
 382         fp->block_offset = 0;
 383     }
 384     fp->block_address = block_address;
 385     fp->block_length = count;
 386     return 0;
 387 }
 388
 389 int
 390 bgzf_read(BGZF* fp, void* data, int length)
 391 {
 392     if (length <= 0) {
 393         return 0;
 394     }
 395     if (fp->open_mode != 'r') {
 396         report_error(fp, "file not open for reading");
 397         return -1;
 398     }
 399
 400     int bytes_read = 0;
 401     byte* output = data;
 402     while (bytes_read < length) {
 403         int available = fp->block_length - fp->block_offset;
 404         if (available <= 0) {
 405             if (read_block(fp) != 0) {
 406                 return -1;
 407             }
 408             available = fp->block_length - fp->block_offset;
 409             if (available <= 0) {
 410                 break;
 411             }
 412         }
 413         int copy_length = min(length-bytes_read, available);
 414         byte* buffer = fp->uncompressed_block;
 415         memcpy(output, buffer + fp->block_offset, copy_length);
 416         fp->block_offset += copy_length;
 417         output += copy_length;
 418         bytes_read += copy_length;
 419     }
 420     if (fp->block_offset == fp->block_length) {
 421 #ifdef _USE_KNETFILE
 422         fp->block_address = knet_tell(fp->x.fpr);
 423 #else
 424         fp->block_address = ftello(fp->file);
 425 #endif
 426         fp->block_offset = 0;
 427         fp->block_length = 0;
 428     }
 429     return bytes_read;
 430 }
 431
 432 static
 433 int
 434 flush_block(BGZF* fp)
 435 {
 436     while (fp->block_offset > 0) {
 437         int block_length = deflate_block(fp, fp->block_offset);
 438         if (block_length < 0) {
 439             return -1;
 440         }
 441 #ifdef _USE_KNETFILE
 442         int count = fwrite(fp->compressed_block, 1, block_length, fp->x.fpw);
 443 #else
 444         int count = fwrite(fp->compressed_block, 1, block_length, fp->file);
 445 #endif
 446         if (count != block_length) {
 447             report_error(fp, "write failed");
 448             return -1;
 449         }
 450         fp->block_address += block_length;
 451     }
 452     return 0;
 453 }
 454
 455 int
 456 bgzf_write(BGZF* fp, const void* data, int length)
 457 {
 458     if (fp->open_mode != 'w') {
 459         report_error(fp, "file not open for writing");
 460         return -1;
 461     }
 462
 463     if (fp->uncompressed_block == NULL) {
 464         fp->uncompressed_block = malloc(fp->uncompressed_block_size);
 465     }
 466
 467     const byte* input = data;
 468     int block_length = fp->uncompressed_block_size;
 469     int bytes_written = 0;
 470     while (bytes_written < length) {
 471         int copy_length = min(block_length - fp->block_offset, length - bytes_written);
 472         byte* buffer = fp->uncompressed_block;
 473         memcpy(buffer + fp->block_offset, input, copy_length);
 474         fp->block_offset += copy_length;
 475         input += copy_length;
 476         bytes_written += copy_length;
 477         if (fp->block_offset == block_length) {
 478             if (flush_block(fp) != 0) {
 479                 break;
 480             }
 481         }
 482     }
 483     return bytes_written;
 484 }
 485
 486 int
 487 bgzf_close(BGZF* fp)
 488 {
 489     if (fp->open_mode == 'w') {
 490         if (flush_block(fp) != 0) {
 491             return -1;
 492         }
 493 #ifdef _USE_KNETFILE
 494         if (fflush(fp->x.fpw) != 0) {
 495 #else
 496         if (fflush(fp->file) != 0) {
 497 #endif
 498             report_error(fp, "flush failed");
 499             return -1;
 500         }
 501     }
 502     if (fp->owned_file) {
 503 #ifdef _USE_KNETFILE
 504                 int ret;
 505                 if (fp->open_mode == 'w') ret = fclose(fp->x.fpw);
 506                 else ret = knet_close(fp->x.fpr);
 507         if (ret != 0) return -1;
 508 #else
 509         if (fclose(fp->file) != 0) {
 510             return -1;
 511         }
 512 #endif
 513     }
 514     free(fp->uncompressed_block);
 515     free(fp->compressed_block);
 516     free(fp);
 517     return 0;
 518 }
 519
 520 int64_t
 521 bgzf_tell(BGZF* fp)
 522 {
 523     return ((fp->block_address << 16) | (fp->block_offset & 0xFFFF));
 524 }
 525
 526 int64_t
 527 bgzf_seek(BGZF* fp, int64_t pos, int where)
 528 {
 529     if (fp->open_mode != 'r') {
 530         report_error(fp, "file not open for read");
 531         return -1;
 532     }
 533     if (where != SEEK_SET) {
 534         report_error(fp, "unimplemented seek option");
 535         return -1;
 536     }
 537     int block_offset = pos & 0xFFFF;
 538     int64_t block_address = (pos >> 16) & 0xFFFFFFFFFFFFLL;
 539 #ifdef _USE_KNETFILE
 540     if (knet_seek(fp->x.fpr, block_address, SEEK_SET) != 0) {
 541 #else
 542     if (fseeko(fp->file, block_address, SEEK_SET) != 0) {
 543 #endif
 544         report_error(fp, "seek failed");
 545         return -1;
 546     }
 547     fp->block_length = 0;  // indicates current block is not loaded
 548     fp->block_address = block_address;
 549     fp->block_offset = block_offset;
 550     return 0;
 551 }
 552