tabix/bgzf.c.pysam.c

   1 #include "pysam.h"
   2
   3 /* The MIT License
   4
   5    Copyright (c) 2008 Broad Institute / Massachusetts Institute of Technology
   6
   7    Permission is hereby granted, free of charge, to any person obtaining a copy
   8    of this software and associated documentation files (the "Software"), to deal
   9    in the Software without restriction, including without limitation the rights
  10    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  11    copies of the Software, and to permit persons to whom the Software is
  12    furnished to do so, subject to the following conditions:
  13
  14    The above copyright notice and this permission notice shall be included in
  15    all copies or substantial portions of the Software.
  16
  17    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  18    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  20    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  21    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  22    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  23    THE SOFTWARE.
  24 */
  25
  26 /*
  27   2009-06-29 by lh3: cache recent uncompressed blocks.
  28   2009-06-25 by lh3: optionally use my knetfile library to access file on a FTP.
  29   2009-06-12 by lh3: support a mode string like "wu" where 'u' for uncompressed output */
  30
  31 #include <stdio.h>
  32 #include <stdlib.h>
  33 #include <string.h>
  34 #include <unistd.h>
  35 #include <fcntl.h>
  36 #include <sys/types.h>
  37 #include <sys/stat.h>
  38 #include "bgzf.h"
  39
  40 #include "khash.h"
  41 typedef struct {
  42         int size;
  43         uint8_t *block;
  44         int64_t end_offset;
  45 } cache_t;
  46 KHASH_MAP_INIT_INT64(cache, cache_t)
  47
  48 #if defined(_WIN32) || defined(_MSC_VER)
  49 #define ftello(fp) ftell(fp)
  50 #define fseeko(fp, offset, whence) fseek(fp, offset, whence)
  51 #else
  52 extern off_t ftello(FILE *stream);
  53 extern int fseeko(FILE *stream, off_t offset, int whence);
  54 #endif
  55
  56 typedef int8_t bgzf_byte_t;
  57
  58 static const int DEFAULT_BLOCK_SIZE = 64 * 1024;
  59 static const int MAX_BLOCK_SIZE = 64 * 1024;
  60
  61 static const int BLOCK_HEADER_LENGTH = 18;
  62 static const int BLOCK_FOOTER_LENGTH = 8;
  63
  64 static const int GZIP_ID1 = 31;
  65 static const int GZIP_ID2 = 139;
  66 static const int CM_DEFLATE = 8;
  67 static const int FLG_FEXTRA = 4;
  68 static const int OS_UNKNOWN = 255;
  69 static const int BGZF_ID1 = 66; // 'B'
  70 static const int BGZF_ID2 = 67; // 'C'
  71 static const int BGZF_LEN = 2;
  72 static const int BGZF_XLEN = 6; // BGZF_LEN+4
  73
  74 static const int GZIP_WINDOW_BITS = -15; // no zlib header
  75 static const int Z_DEFAULT_MEM_LEVEL = 8;
  76
  77
  78 inline
  79 void
  80 packInt16(uint8_t* buffer, uint16_t value)
  81 {
  82     buffer[0] = value;
  83     buffer[1] = value >> 8;
  84 }
  85
  86 inline
  87 int
  88 unpackInt16(const uint8_t* buffer)
  89 {
  90     return (buffer[0] | (buffer[1] << 8));
  91 }
  92
  93 inline
  94 void
  95 packInt32(uint8_t* buffer, uint32_t value)
  96 {
  97     buffer[0] = value;
  98     buffer[1] = value >> 8;
  99     buffer[2] = value >> 16;
 100     buffer[3] = value >> 24;
 101 }
 102
 103 static inline
 104 int
 105 bgzf_min(int x, int y)
 106 {
 107     return (x < y) ? x : y;
 108 }
 109
 110 static
 111 void
 112 report_error(BGZF* fp, const char* message) {
 113     fp->error = message;
 114 }
 115
 116 int bgzf_check_bgzf(const char *fn)
 117 {
 118     BGZF *fp;
 119     uint8_t buf[10],magic[10]="\037\213\010\4\0\0\0\0\0\377";
 120     int n;
 121
 122     if ((fp = bgzf_open(fn, "r")) == 0)
 123     {
 124         fprintf(pysamerr, "[bgzf_check_bgzf] failed to open the file: %s\n",fn);
 125         return -1;
 126     }
 127
 128 #ifdef _USE_KNETFILE
 129     n = knet_read(fp->x.fpr, buf, 10);
 130 #else
 131     n = fread(buf, 1, 10, fp->file);
 132 #endif
 133     bgzf_close(fp);
 134
 135     if ( n!=10 )
 136         return -1;
 137
 138     if ( !memcmp(magic, buf, 10) ) return 1;
 139     return 0;
 140 }
 141
 142 static BGZF *bgzf_read_init()
 143 {
 144         BGZF *fp;
 145         fp = calloc(1, sizeof(BGZF));
 146     fp->uncompressed_block_size = MAX_BLOCK_SIZE;
 147     fp->uncompressed_block = malloc(MAX_BLOCK_SIZE);
 148     fp->compressed_block_size = MAX_BLOCK_SIZE;
 149     fp->compressed_block = malloc(MAX_BLOCK_SIZE);
 150         fp->cache_size = 0;
 151         fp->cache = kh_init(cache);
 152         return fp;
 153 }
 154
 155 static
 156 BGZF*
 157 open_read(int fd)
 158 {
 159 #ifdef _USE_KNETFILE
 160     knetFile *file = knet_dopen(fd, "r");
 161 #else
 162     FILE* file = fdopen(fd, "r");
 163 #endif
 164     BGZF* fp;
 165         if (file == 0) return 0;
 166         fp = bgzf_read_init();
 167     fp->file_descriptor = fd;
 168     fp->open_mode = 'r';
 169 #ifdef _USE_KNETFILE
 170     fp->x.fpr = file;
 171 #else
 172     fp->file = file;
 173 #endif
 174     return fp;
 175 }
 176
 177 static
 178 BGZF*
 179 open_write(int fd, int compress_level) // compress_level==-1 for the default level
 180 {
 181     FILE* file = fdopen(fd, "w");
 182     BGZF* fp;
 183         if (file == 0) return 0;
 184         fp = malloc(sizeof(BGZF));
 185     fp->file_descriptor = fd;
 186     fp->open_mode = 'w';
 187     fp->owned_file = 0;
 188         fp->compress_level = compress_level < 0? Z_DEFAULT_COMPRESSION : compress_level; // Z_DEFAULT_COMPRESSION==-1
 189         if (fp->compress_level > 9) fp->compress_level = Z_DEFAULT_COMPRESSION;
 190 #ifdef _USE_KNETFILE
 191     fp->x.fpw = file;
 192 #else
 193     fp->file = file;
 194 #endif
 195     fp->uncompressed_block_size = DEFAULT_BLOCK_SIZE;
 196     fp->uncompressed_block = NULL;
 197     fp->compressed_block_size = MAX_BLOCK_SIZE;
 198     fp->compressed_block = malloc(MAX_BLOCK_SIZE);
 199     fp->block_address = 0;
 200     fp->block_offset = 0;
 201     fp->block_length = 0;
 202     fp->error = NULL;
 203     return fp;
 204 }
 205
 206 BGZF*
 207 bgzf_open(const char* __restrict path, const char* __restrict mode)
 208 {
 209     BGZF* fp = NULL;
 210     if (strchr(mode, 'r') || strchr(mode, 'R')) { /* The reading mode is preferred. */
 211 #ifdef _USE_KNETFILE
 212                 knetFile *file = knet_open(path, mode);
 213                 if (file == 0) return 0;
 214                 fp = bgzf_read_init();
 215                 fp->file_descriptor = -1;
 216                 fp->open_mode = 'r';
 217                 fp->x.fpr = file;
 218 #else
 219                 int fd, oflag = O_RDONLY;
 220 #ifdef _WIN32
 221                 oflag |= O_BINARY;
 222 #endif
 223                 fd = open(path, oflag);
 224                 if (fd == -1) return 0;
 225         fp = open_read(fd);
 226 #endif
 227     } else if (strchr(mode, 'w') || strchr(mode, 'W')) {
 228                 int fd, compress_level = -1, oflag = O_WRONLY | O_CREAT | O_TRUNC;
 229 #ifdef _WIN32
 230                 oflag |= O_BINARY;
 231 #endif
 232                 fd = open(path, oflag, 0666);
 233                 if (fd == -1) return 0;
 234                 { // set compress_level
 235                         int i;
 236                         for (i = 0; mode[i]; ++i)
 237                                 if (mode[i] >= '0' && mode[i] <= '9') break;
 238                         if (mode[i]) compress_level = (int)mode[i] - '0';
 239                         if (strchr(mode, 'u')) compress_level = 0;
 240                 }
 241         fp = open_write(fd, compress_level);
 242     }
 243     if (fp != NULL) fp->owned_file = 1;
 244     return fp;
 245 }
 246
 247 BGZF*
 248 bgzf_fdopen(int fd, const char * __restrict mode)
 249 {
 250         if (fd == -1) return 0;
 251     if (mode[0] == 'r' || mode[0] == 'R') {
 252         return open_read(fd);
 253     } else if (mode[0] == 'w' || mode[0] == 'W') {
 254                 int i, compress_level = -1;
 255                 for (i = 0; mode[i]; ++i)
 256                         if (mode[i] >= '0' && mode[i] <= '9') break;
 257                 if (mode[i]) compress_level = (int)mode[i] - '0';
 258                 if (strchr(mode, 'u')) compress_level = 0;
 259         return open_write(fd, compress_level);
 260     } else {
 261         return NULL;
 262     }
 263 }
 264
 265 static
 266 int
 267 deflate_block(BGZF* fp, int block_length)
 268 {
 269     // Deflate the block in fp->uncompressed_block into fp->compressed_block.
 270     // Also adds an extra field that stores the compressed block length.
 271
 272     bgzf_byte_t* buffer = fp->compressed_block;
 273     int buffer_size = fp->compressed_block_size;
 274
 275     // Init gzip header
 276     buffer[0] = GZIP_ID1;
 277     buffer[1] = GZIP_ID2;
 278     buffer[2] = CM_DEFLATE;
 279     buffer[3] = FLG_FEXTRA;
 280     buffer[4] = 0; // mtime
 281     buffer[5] = 0;
 282     buffer[6] = 0;
 283     buffer[7] = 0;
 284     buffer[8] = 0;
 285     buffer[9] = OS_UNKNOWN;
 286     buffer[10] = BGZF_XLEN;
 287     buffer[11] = 0;
 288     buffer[12] = BGZF_ID1;
 289     buffer[13] = BGZF_ID2;
 290     buffer[14] = BGZF_LEN;
 291     buffer[15] = 0;
 292     buffer[16] = 0; // placeholder for block length
 293     buffer[17] = 0;
 294
 295     // loop to retry for blocks that do not compress enough
 296     int input_length = block_length;
 297     int compressed_length = 0;
 298     while (1) {
 299         z_stream zs;
 300         zs.zalloc = NULL;
 301         zs.zfree = NULL;
 302         zs.next_in = fp->uncompressed_block;
 303         zs.avail_in = input_length;
 304         zs.next_out = (void*)&buffer[BLOCK_HEADER_LENGTH];
 305         zs.avail_out = buffer_size - BLOCK_HEADER_LENGTH - BLOCK_FOOTER_LENGTH;
 306
 307         int status = deflateInit2(&zs, fp->compress_level, Z_DEFLATED,
 308                                   GZIP_WINDOW_BITS, Z_DEFAULT_MEM_LEVEL, Z_DEFAULT_STRATEGY);
 309         if (status != Z_OK) {
 310             report_error(fp, "deflate init failed");
 311             return -1;
 312         }
 313         status = deflate(&zs, Z_FINISH);
 314         if (status != Z_STREAM_END) {
 315             deflateEnd(&zs);
 316             if (status == Z_OK) {
 317                 // Not enough space in buffer.
 318                 // Can happen in the rare case the input doesn't compress enough.
 319                 // Reduce the amount of input until it fits.
 320                 input_length -= 1024;
 321                 if (input_length <= 0) {
 322                     // should never happen
 323                     report_error(fp, "input reduction failed");
 324                     return -1;
 325                 }
 326                 continue;
 327             }
 328             report_error(fp, "deflate failed");
 329             return -1;
 330         }
 331         status = deflateEnd(&zs);
 332         if (status != Z_OK) {
 333             report_error(fp, "deflate end failed");
 334             return -1;
 335         }
 336         compressed_length = zs.total_out;
 337         compressed_length += BLOCK_HEADER_LENGTH + BLOCK_FOOTER_LENGTH;
 338         if (compressed_length > MAX_BLOCK_SIZE) {
 339             // should never happen
 340             report_error(fp, "deflate overflow");
 341             return -1;
 342         }
 343         break;
 344     }
 345
 346     packInt16((uint8_t*)&buffer[16], compressed_length-1);
 347     uint32_t crc = crc32(0L, NULL, 0L);
 348     crc = crc32(crc, fp->uncompressed_block, input_length);
 349     packInt32((uint8_t*)&buffer[compressed_length-8], crc);
 350     packInt32((uint8_t*)&buffer[compressed_length-4], input_length);
 351
 352     int remaining = block_length - input_length;
 353     if (remaining > 0) {
 354         if (remaining > input_length) {
 355             // should never happen (check so we can use memcpy)
 356             report_error(fp, "remainder too large");
 357             return -1;
 358         }
 359         memcpy(fp->uncompressed_block,
 360                fp->uncompressed_block + input_length,
 361                remaining);
 362     }
 363     fp->block_offset = remaining;
 364     return compressed_length;
 365 }
 366
 367 static
 368 int
 369 inflate_block(BGZF* fp, int block_length)
 370 {
 371     // Inflate the block in fp->compressed_block into fp->uncompressed_block
 372
 373     z_stream zs;
 374         int status;
 375     zs.zalloc = NULL;
 376     zs.zfree = NULL;
 377     zs.next_in = fp->compressed_block + 18;
 378     zs.avail_in = block_length - 16;
 379     zs.next_out = fp->uncompressed_block;
 380     zs.avail_out = fp->uncompressed_block_size;
 381
 382     status = inflateInit2(&zs, GZIP_WINDOW_BITS);
 383     if (status != Z_OK) {
 384         report_error(fp, "inflate init failed");
 385         return -1;
 386     }
 387     status = inflate(&zs, Z_FINISH);
 388     if (status != Z_STREAM_END) {
 389         inflateEnd(&zs);
 390         report_error(fp, "inflate failed");
 391         return -1;
 392     }
 393     status = inflateEnd(&zs);
 394     if (status != Z_OK) {
 395         report_error(fp, "inflate failed");
 396         return -1;
 397     }
 398     return zs.total_out;
 399 }
 400
 401 static
 402 int
 403 check_header(const bgzf_byte_t* header)
 404 {
 405     return (header[0] == GZIP_ID1 &&
 406             header[1] == (bgzf_byte_t) GZIP_ID2 &&
 407             header[2] == Z_DEFLATED &&
 408             (header[3] & FLG_FEXTRA) != 0 &&
 409             unpackInt16((uint8_t*)&header[10]) == BGZF_XLEN &&
 410             header[12] == BGZF_ID1 &&
 411             header[13] == BGZF_ID2 &&
 412             unpackInt16((uint8_t*)&header[14]) == BGZF_LEN);
 413 }
 414
 415 static void free_cache(BGZF *fp)
 416 {
 417         khint_t k;
 418         khash_t(cache) *h = (khash_t(cache)*)fp->cache;
 419         if (fp->open_mode != 'r') return;
 420         for (k = kh_begin(h); k < kh_end(h); ++k)
 421                 if (kh_exist(h, k)) free(kh_val(h, k).block);
 422         kh_destroy(cache, h);
 423 }
 424
 425 static int load_block_from_cache(BGZF *fp, int64_t block_address)
 426 {
 427         khint_t k;
 428         cache_t *p;
 429         khash_t(cache) *h = (khash_t(cache)*)fp->cache;
 430         k = kh_get(cache, h, block_address);
 431         if (k == kh_end(h)) return 0;
 432         p = &kh_val(h, k);
 433         if (fp->block_length != 0) fp->block_offset = 0;
 434         fp->block_address = block_address;
 435         fp->block_length = p->size;
 436         memcpy(fp->uncompressed_block, p->block, MAX_BLOCK_SIZE);
 437 #ifdef _USE_KNETFILE
 438         knet_seek(fp->x.fpr, p->end_offset, SEEK_SET);
 439 #else
 440         fseeko(fp->file, p->end_offset, SEEK_SET);
 441 #endif
 442         return p->size;
 443 }
 444
 445 static void cache_block(BGZF *fp, int size)
 446 {
 447         int ret;
 448         khint_t k;
 449         cache_t *p;
 450         khash_t(cache) *h = (khash_t(cache)*)fp->cache;
 451         if (MAX_BLOCK_SIZE >= fp->cache_size) return;
 452         if ((kh_size(h) + 1) * MAX_BLOCK_SIZE > fp->cache_size) {
 453                 /* A better way would be to remove the oldest block in the
 454                  * cache, but here we remove a random one for simplicity. This
 455                  * should not have a big impact on performance. */
 456                 for (k = kh_begin(h); k < kh_end(h); ++k)
 457                         if (kh_exist(h, k)) break;
 458                 if (k < kh_end(h)) {
 459                         free(kh_val(h, k).block);
 460                         kh_del(cache, h, k);
 461                 }
 462         }
 463         k = kh_put(cache, h, fp->block_address, &ret);
 464         if (ret == 0) return; // if this happens, a bug!
 465         p = &kh_val(h, k);
 466         p->size = fp->block_length;
 467         p->end_offset = fp->block_address + size;
 468         p->block = malloc(MAX_BLOCK_SIZE);
 469         memcpy(kh_val(h, k).block, fp->uncompressed_block, MAX_BLOCK_SIZE);
 470 }
 471
 472 int
 473 bgzf_read_block(BGZF* fp)
 474 {
 475     bgzf_byte_t header[BLOCK_HEADER_LENGTH];
 476         int count, size = 0, block_length, remaining;
 477 #ifdef _USE_KNETFILE
 478     int64_t block_address = knet_tell(fp->x.fpr);
 479         if (load_block_from_cache(fp, block_address)) return 0;
 480     count = knet_read(fp->x.fpr, header, sizeof(header));
 481 #else
 482     int64_t block_address = ftello(fp->file);
 483         if (load_block_from_cache(fp, block_address)) return 0;
 484     count = fread(header, 1, sizeof(header), fp->file);
 485 #endif
 486     if (count == 0) {
 487         fp->block_length = 0;
 488         return 0;
 489     }
 490         size = count;
 491     if (count != sizeof(header)) {
 492         report_error(fp, "read failed");
 493         return -1;
 494     }
 495     if (!check_header(header)) {
 496         report_error(fp, "invalid block header");
 497         return -1;
 498     }
 499     block_length = unpackInt16((uint8_t*)&header[16]) + 1;
 500     bgzf_byte_t* compressed_block = (bgzf_byte_t*) fp->compressed_block;
 501     memcpy(compressed_block, header, BLOCK_HEADER_LENGTH);
 502     remaining = block_length - BLOCK_HEADER_LENGTH;
 503 #ifdef _USE_KNETFILE
 504     count = knet_read(fp->x.fpr, &compressed_block[BLOCK_HEADER_LENGTH], remaining);
 505 #else
 506     count = fread(&compressed_block[BLOCK_HEADER_LENGTH], 1, remaining, fp->file);
 507 #endif
 508     if (count != remaining) {
 509         report_error(fp, "read failed");
 510         return -1;
 511     }
 512         size += count;
 513     count = inflate_block(fp, block_length);
 514     if (count < 0) return -1;
 515     if (fp->block_length != 0) {
 516         // Do not reset offset if this read follows a seek.
 517         fp->block_offset = 0;
 518     }
 519     fp->block_address = block_address;
 520     fp->block_length = count;
 521         cache_block(fp, size);
 522     return 0;
 523 }
 524
 525 int
 526 bgzf_read(BGZF* fp, void* data, int length)
 527 {
 528     if (length <= 0) {
 529         return 0;
 530     }
 531     if (fp->open_mode != 'r') {
 532         report_error(fp, "file not open for reading");
 533         return -1;
 534     }
 535
 536     int bytes_read = 0;
 537     bgzf_byte_t* output = data;
 538     while (bytes_read < length) {
 539         int copy_length, available = fp->block_length - fp->block_offset;
 540                 bgzf_byte_t *buffer;
 541         if (available <= 0) {
 542             if (bgzf_read_block(fp) != 0) {
 543                 return -1;
 544             }
 545             available = fp->block_length - fp->block_offset;
 546             if (available <= 0) {
 547                 break;
 548             }
 549         }
 550         copy_length = bgzf_min(length-bytes_read, available);
 551         buffer = fp->uncompressed_block;
 552         memcpy(output, buffer + fp->block_offset, copy_length);
 553         fp->block_offset += copy_length;
 554         output += copy_length;
 555         bytes_read += copy_length;
 556     }
 557     if (fp->block_offset == fp->block_length) {
 558 #ifdef _USE_KNETFILE
 559         fp->block_address = knet_tell(fp->x.fpr);
 560 #else
 561         fp->block_address = ftello(fp->file);
 562 #endif
 563         fp->block_offset = 0;
 564         fp->block_length = 0;
 565     }
 566     return bytes_read;
 567 }
 568
 569 int bgzf_flush(BGZF* fp)
 570 {
 571     while (fp->block_offset > 0) {
 572         int count, block_length;
 573                 block_length = deflate_block(fp, fp->block_offset);
 574         if (block_length < 0) return -1;
 575 #ifdef _USE_KNETFILE
 576         count = fwrite(fp->compressed_block, 1, block_length, fp->x.fpw);
 577 #else
 578         count = fwrite(fp->compressed_block, 1, block_length, fp->file);
 579 #endif
 580         if (count != block_length) {
 581             report_error(fp, "write failed");
 582             return -1;
 583         }
 584         fp->block_address += block_length;
 585     }
 586     return 0;
 587 }
 588
 589 int bgzf_flush_try(BGZF *fp, int size)
 590 {
 591         if (fp->block_offset + size > fp->uncompressed_block_size)
 592                 return bgzf_flush(fp);
 593         return -1;
 594 }
 595
 596 int bgzf_write(BGZF* fp, const void* data, int length)
 597 {
 598         const bgzf_byte_t *input = data;
 599         int block_length, bytes_written;
 600     if (fp->open_mode != 'w') {
 601         report_error(fp, "file not open for writing");
 602         return -1;
 603     }
 604
 605     if (fp->uncompressed_block == NULL)
 606         fp->uncompressed_block = malloc(fp->uncompressed_block_size);
 607
 608     input = data;
 609     block_length = fp->uncompressed_block_size;
 610     bytes_written = 0;
 611     while (bytes_written < length) {
 612         int copy_length = bgzf_min(block_length - fp->block_offset, length - bytes_written);
 613         bgzf_byte_t* buffer = fp->uncompressed_block;
 614         memcpy(buffer + fp->block_offset, input, copy_length);
 615         fp->block_offset += copy_length;
 616         input += copy_length;
 617         bytes_written += copy_length;
 618         if (fp->block_offset == block_length) {
 619             if (bgzf_flush(fp) != 0) {
 620                 break;
 621             }
 622         }
 623     }
 624     return bytes_written;
 625 }
 626
 627 int bgzf_close(BGZF* fp)
 628 {
 629     if (fp->open_mode == 'w') {
 630         if (bgzf_flush(fp) != 0) return -1;
 631                 { // add an empty block
 632                         int count, block_length = deflate_block(fp, 0);
 633 #ifdef _USE_KNETFILE
 634                         count = fwrite(fp->compressed_block, 1, block_length, fp->x.fpw);
 635 #else
 636                         count = fwrite(fp->compressed_block, 1, block_length, fp->file);
 637 #endif
 638                 }
 639 #ifdef _USE_KNETFILE
 640         if (fflush(fp->x.fpw) != 0) {
 641 #else
 642         if (fflush(fp->file) != 0) {
 643 #endif
 644             report_error(fp, "flush failed");
 645             return -1;
 646         }
 647     }
 648     if (fp->owned_file) {
 649 #ifdef _USE_KNETFILE
 650                 int ret;
 651                 if (fp->open_mode == 'w') ret = fclose(fp->x.fpw);
 652                 else ret = knet_close(fp->x.fpr);
 653         if (ret != 0) return -1;
 654 #else
 655         if (fclose(fp->file) != 0) return -1;
 656 #endif
 657     }
 658     free(fp->uncompressed_block);
 659     free(fp->compressed_block);
 660         free_cache(fp);
 661     free(fp);
 662     return 0;
 663 }
 664
 665 void bgzf_set_cache_size(BGZF *fp, int cache_size)
 666 {
 667         if (fp) fp->cache_size = cache_size;
 668 }
 669
 670 int bgzf_check_EOF(BGZF *fp)
 671 {
 672         static uint8_t magic[28] = "\037\213\010\4\0\0\0\0\0\377\6\0\102\103\2\0\033\0\3\0\0\0\0\0\0\0\0\0";
 673         uint8_t buf[28];
 674         off_t offset;
 675 #ifdef _USE_KNETFILE
 676         offset = knet_tell(fp->x.fpr);
 677         if (knet_seek(fp->x.fpr, -28, SEEK_END) != 0) return -1;
 678         knet_read(fp->x.fpr, buf, 28);
 679         knet_seek(fp->x.fpr, offset, SEEK_SET);
 680 #else
 681         offset = ftello(fp->file);
 682         if (fseeko(fp->file, -28, SEEK_END) != 0) return -1;
 683         fread(buf, 1, 28, fp->file);
 684         fseeko(fp->file, offset, SEEK_SET);
 685 #endif
 686         return (memcmp(magic, buf, 28) == 0)? 1 : 0;
 687 }
 688
 689 int64_t bgzf_seek(BGZF* fp, int64_t pos, int where)
 690 {
 691         int block_offset;
 692         int64_t block_address;
 693
 694     if (fp->open_mode != 'r') {
 695         report_error(fp, "file not open for read");
 696         return -1;
 697     }
 698     if (where != SEEK_SET) {
 699         report_error(fp, "unimplemented seek option");
 700         return -1;
 701     }
 702     block_offset = pos & 0xFFFF;
 703     block_address = (pos >> 16) & 0xFFFFFFFFFFFFLL;
 704 #ifdef _USE_KNETFILE
 705     if (knet_seek(fp->x.fpr, block_address, SEEK_SET) != 0) {
 706 #else
 707     if (fseeko(fp->file, block_address, SEEK_SET) != 0) {
 708 #endif
 709         report_error(fp, "seek failed");
 710         return -1;
 711     }
 712     fp->block_length = 0;  // indicates current block is not loaded
 713     fp->block_address = block_address;
 714     fp->block_offset = block_offset;
 715     return 0;
 716 }