text: /* -*- mode: javascript; c-basic-offset: 4; indent-tabs-mode: nil -*- */
  
  // 
  // Dalliance Genome Explorer
  // (c) Thomas Down 2006-2011
  //
  // bam.js: indexed binary alignments
  //
  
  "use strict";
  
  var BAM_MAGIC = 21840194;
  var BAI_MAGIC = 21578050;
  
  function BamFile() {
  }
  
  function Vob(b, o) {
      this.block = b;
      this.offset = o;
  }
  
  Vob.prototype.toString = function() {
      return '' + this.block + ':' + this.offset;
  }
  
  function Chunk(minv, maxv) {
      this.minv = minv; this.maxv = maxv;
  }
  
  function makeBam(data, bai, callback) {
      var bam = new BamFile();
      bam.data = data;
      bam.bai = bai;
  
      bam.bai.fetch(function(header) {   // Do we really need to fetch the whole thing? :-(
          if (!header) {
              return callback(null, "Couldn't access BAI");
          }
  
          var uncba = new Uint8Array(header);
          var baiMagic = readInt(uncba, 0);
          if (baiMagic != BAI_MAGIC) {
              return callback(null, 'Not a BAI file, magic=0x' + baiMagic.toString(16));
          }
  
          var nref = readInt(uncba, 4);
  
          bam.indices = [];
  
          var p = 8;
          var minBlockIndex = 1000000000;
          for (var ref = 0; ref < nref; ++ref) {
              var blockStart = p;
              var nbin = readInt(uncba, p); p += 4;
              for (var b = 0; b < nbin; ++b) {
                  var bin = readInt(uncba, p);
                  var nchnk = readInt(uncba, p+4);
                  p += 8 + (nchnk * 16);
              }
              var nintv = readInt(uncba, p); p += 4;
              
              var q = p;
              for (var i = 0; i < nintv; ++i) {
                  var v = readVob(uncba, q); q += 8;
                  if (v) {
                      var bi = v.block;
                      if (v.offset > 0)
                          bi += 65536;
  
                      if (bi < minBlockIndex)
                          minBlockIndex = bi;
                      break;
                  }
              }
              p += (nintv * 8);
  
  
              if (nbin > 0) {
                  bam.indices[ref] = new Uint8Array(header, blockStart, p - blockStart);
              }                     
          }
  
          bam.data.slice(0, minBlockIndex).fetch(function(r) {
              if (!r) {
                  return callback(null, "Couldn't access BAM");
              }
              
              var unc = unbgzf(r, r.byteLength);
              var uncba = new Uint8Array(unc);
  
              var magic = readInt(uncba, 0);
              if (magic != BAM_MAGIC) {
                  return callback(null, "Not a BAM file, magic=0x" + magic.toString(16));
              }
              var headLen = readInt(uncba, 4);
              var header = '';
              for (var i = 0; i < headLen; ++i) {
                  header += String.fromCharCode(uncba[i + 8]);
              }
  
              var nRef = readInt(uncba, headLen + 8);
              var p = headLen + 12;
  
              bam.chrToIndex = {};
              bam.indexToChr = [];
              for (var i = 0; i < nRef; ++i) {
                  var lName = readInt(uncba, p);
                  var name = '';
                  for (var j = 0; j < lName-1; ++j) {
                      name += String.fromCharCode(uncba[p + 4 + j]);
                  }
                  var lRef = readInt(uncba, p + lName + 4);
                  bam.chrToIndex[name] = i;
                  if (name.indexOf('chr') == 0) {
                      bam.chrToIndex[name.substring(3)] = i;
                  } else {
                      bam.chrToIndex['chr' + name] = i;
                  }
                  bam.indexToChr.push(name);
  
                  p = p + 8 + lName;
              }
  
              if (bam.indices) {
                  return callback(bam);
              }
          });
      });
  }
  
  
  
  BamFile.prototype.blocksForRange = function(refId, min, max) {
      var index = this.indices[refId];
      if (!index) {
          return [];
      }
  
      var intBinsL = reg2bins(min, max);
      var intBins = [];
      for (var i = 0; i < intBinsL.length; ++i) {
          intBins[intBinsL[i]] = true;
      }
      var leafChunks = [], otherChunks = [];
  
      var nbin = readInt(index, 0);
      var p = 4;
      for (var b = 0; b < nbin; ++b) {
          var bin = readInt(index, p);
          var nchnk = readInt(index, p+4);
  //        dlog('bin=' + bin + '; nchnk=' + nchnk);
          p += 8;
          if (intBins[bin]) {
              for (var c = 0; c < nchnk; ++c) {
                  var cs = readVob(index, p);
                  var ce = readVob(index, p + 8);
                  (bin < 4681 ? otherChunks : leafChunks).push(new Chunk(cs, ce));
                  p += 16;
              }
          } else {
              p +=  (nchnk * 16);
          }
      }
  //    dlog('leafChunks = ' + miniJSONify(leafChunks));
  //    dlog('otherChunks = ' + miniJSONify(otherChunks));
  
      var nintv = readInt(index, p);
      var lowest = null;
      var minLin = Math.min(min>>14, nintv - 1), maxLin = Math.min(max>>14, nintv - 1);
      for (var i = minLin; i <= maxLin; ++i) {
          var lb =  readVob(index, p + 4 + (i * 8));
          if (!lb) {
              continue;
          }
          if (!lowest || lb.block < lowest.block || lb.offset < lowest.offset) {
              lowest = lb;
          }
      }
      // dlog('Lowest LB = ' + lowest);
      
      var prunedOtherChunks = [];
      if (lowest != null) {
          for (var i = 0; i < otherChunks.length; ++i) {
              var chnk = otherChunks[i];
              if (chnk.maxv.block >= lowest.block && chnk.maxv.offset >= lowest.offset) {
                  prunedOtherChunks.push(chnk);
              }
          }
      }
      // dlog('prunedOtherChunks = ' + miniJSONify(prunedOtherChunks));
      otherChunks = prunedOtherChunks;
  
      var intChunks = [];
      for (var i = 0; i < otherChunks.length; ++i) {
          intChunks.push(otherChunks[i]);
      }
      for (var i = 0; i < leafChunks.length; ++i) {
          intChunks.push(leafChunks[i]);
      }
  
      intChunks.sort(function(c0, c1) {
          var dif = c0.minv.block - c1.minv.block;
          if (dif != 0) {
              return dif;
          } else {
              return c0.minv.offset - c1.minv.offset;
          }
      });
      var mergedChunks = [];
      if (intChunks.length > 0) {
          var cur = intChunks[0];
          for (var i = 1; i < intChunks.length; ++i) {
              var nc = intChunks[i];
              if (nc.minv.block == cur.maxv.block /* && nc.minv.offset == cur.maxv.offset */) { // no point splitting mid-block
                  cur = new Chunk(cur.minv, nc.maxv);
              } else {
                  mergedChunks.push(cur);
                  cur = nc;
              }
          }
          mergedChunks.push(cur);
      }
      // dlog('mergedChunks = ' + miniJSONify(mergedChunks));
  
      return mergedChunks;
  }
  
  BamFile.prototype.fetch = function(chr, min, max, callback, opts) {
      var thisB = this;
      opts = opts || {};
  
      var chrId = this.chrToIndex[chr];
      var chunks;
      if (chrId === undefined) {
          chunks = [];
      } else {
          chunks = this.blocksForRange(chrId, min, max);
          if (!chunks) {
              callback(null, 'Error in index fetch');
          }
      }
      
      var records = [];
      var index = 0;
      var data;
  
      function tramp() {
          if (index >= chunks.length) {
              return callback(records);
          } else if (!data) {
              // dlog('fetching ' + index);
              var c = chunks[index];
              var fetchMin = c.minv.block;
              var fetchMax = c.maxv.block + (1<<16); // *sigh*
              thisB.data.slice(fetchMin, fetchMax - fetchMin).fetch(function(r) {
                  data = unbgzf(r, c.maxv.block - c.minv.block + 1);
                  return tramp();
              });
          } else {
              var ba = new Uint8Array(data);
              thisB.readBamRecords(ba, chunks[index].minv.offset, records, min, max, chrId, opts);
              data = null;
              ++index;
              return tramp();
          }
      }
      tramp();
  }
  
  var SEQRET_DECODER = ['=', 'A', 'C', 'x', 'G', 'x', 'x', 'x', 'T', 'x', 'x', 'x', 'x', 'x', 'x', 'N'];
  var CIGAR_DECODER = ['M', 'I', 'D', 'N', 'S', 'H', 'P', '=', 'X', '?', '?', '?', '?', '?', '?', '?'];
  
  function BamRecord() {
  }
  
  BamFile.prototype.readBamRecords = function(ba, offset, sink, min, max, chrId, opts) {
      while (true) {
          var blockSize = readInt(ba, offset);
          var blockEnd = offset + blockSize + 4;
          if (blockEnd >= ba.length) {
              return sink;
          }
  
          var record = new BamRecord();
  
          var refID = readInt(ba, offset + 4);
          var pos = readInt(ba, offset + 8);
          
          var bmn = readInt(ba, offset + 12);
          var bin = (bmn & 0xffff0000) >> 16;
          var mq = (bmn & 0xff00) >> 8;
          var nl = bmn & 0xff;
  
          var flag_nc = readInt(ba, offset + 16);
          var flag = (flag_nc & 0xffff0000) >> 16;
          var nc = flag_nc & 0xffff;
      
          var lseq = readInt(ba, offset + 20);
          
          var nextRef  = readInt(ba, offset + 24);
          var nextPos = readInt(ba, offset + 28);
          
          var tlen = readInt(ba, offset + 32);
      
          record.segment = this.indexToChr[refID];
          record.flag = flag;
          record.pos = pos;
          record.mq = mq;
          if (opts.light)
              record.seqLength = lseq;
  
          if (!opts.light) {
              var readName = '';
              for (var j = 0; j < nl-1; ++j) {
                  readName += String.fromCharCode(ba[offset + 36 + j]);
              }
              record.readName = readName;
          
              var p = offset + 36 + nl;
  
              var cigar = '';
              for (var c = 0; c < nc; ++c) {
                  var cigop = readInt(ba, p);
                  cigar = cigar + (cigop>>4) + CIGAR_DECODER[cigop & 0xf];
                  p += 4;
              }
              record.cigar = cigar;
          
              var seq = '';
              var seqBytes = (lseq + 1) >> 1;
              for (var j = 0; j < seqBytes; ++j) {
                  var sb = ba[p + j];
                  seq += SEQRET_DECODER[(sb & 0xf0) >> 4];
                  seq += SEQRET_DECODER[(sb & 0x0f)];
              }
              p += seqBytes;
              record.seq = seq;
  
              var qseq = '';
              for (var j = 0; j < lseq; ++j) {
                  qseq += String.fromCharCode(ba[p + j] + 33);
              }
              p += lseq;
              record.quals = qseq;
  
              while (p < blockEnd) {
                  var tag = String.fromCharCode(ba[p], ba[p + 1]);
                  var type = String.fromCharCode(ba[p + 2]);
                  var value;
  
                  if (type == 'A') {
                      value = String.fromCharCode(ba[p + 3]);
                      p += 4;
                  } else if (type == 'i' || type == 'I') {
                      value = readInt(ba, p + 3);
                      p += 7;
                  } else if (type == 'c' || type == 'C') {
                      value = ba[p + 3];
                      p += 4;
                  } else if (type == 's' || type == 'S') {
                      value = readShort(ba, p + 3);
                      p += 5;
                  } else if (type == 'f') {
                      value = readFloat(ba, p + 3);
                      p += 7;
                  } else if (type == 'Z' || type == 'H') {
                      p += 3;
                      value = '';
                      for (;;) {
                          var cc = ba[p++];
                          if (cc == 0) {
                              break;
                          } else {
                              value += String.fromCharCode(cc);
                          }
                      }
                  } else if (type == 'B') {
                      var atype = String.fromCharCode(ba[p + 3]);
                      var alen = readInt(ba, p + 4);
                      var elen;
                      var reader;
                      if (atype == 'i' || atype == 'I' || atype == 'f') {
                          elen = 4;
                          if (atype == 'f')
                              reader = readFloat;
                          else
                              reader = readInt;
                      } else if (atype == 's' || atype == 'S') {
                          elen = 2;
                          reader = readShort;
                      } else if (atype == 'c' || atype == 'C') {
                          elen = 1;
                          reader = readByte;
                      } else {
                          throw 'Unknown array type ' + atype;
                      }
  
                      p += 8;
                      value = [];
                      for (var i = 0; i < alen; ++i) {
                          value.push(reader(ba, p));
                          p += elen;
                      }
                  } else {
                      throw 'Unknown type '+ type;
                  }
                  record[tag] = value;
              }
          }
  
          if (!min || record.pos <= max && record.pos + lseq >= min) {
              if (chrId === undefined || refID == chrId) {
                  sink.push(record);
              }
          }
          offset = blockEnd;
      }
  
      // Exits via top of loop.
  };
  
  (function(global) {
      var convertBuffer = new ArrayBuffer(8);
      var ba = new Uint8Array(convertBuffer);
      var fa = new Float32Array(convertBuffer);
  
  
      global.readFloat = function(buf, offset) {
          ba[0] = buf[offset];
          ba[1] = buf[offset+1];
          ba[2] = buf[offset+2];
          ba[3] = buf[offset+3];
          return fa[0];
      };
   }(this));
  
  function readInt64(ba, offset) {
      return (ba[offset + 7] << 24) | (ba[offset + 6] << 16) | (ba[offset + 5] << 8) | (ba[offset + 4]);
  }
  
  function readInt(ba, offset) {
      return (ba[offset + 3] << 24) | (ba[offset + 2] << 16) | (ba[offset + 1] << 8) | (ba[offset]);
  }
  
  function readIntBE(ba, offset) {
      return (ba[offset] << 24) | (ba[offset + 1] << 16) | (ba[offset + 2] << 8) | (ba[offset + 3]);
  }
  
  function readShort(ba, offset) {
      return (ba[offset + 1] << 8) | (ba[offset]);
  }
  
  function readByte(ba, offset) {
      return ba[offset];
  }
  
  function readVob(ba, offset) {
      var block = ((ba[offset+6] & 0xff) * 0x100000000) + ((ba[offset+5] & 0xff) * 0x1000000) + ((ba[offset+4] & 0xff) * 0x10000) + ((ba[offset+3] & 0xff) * 0x100) + ((ba[offset+2] & 0xff));
      var bint = (ba[offset+1] << 8) | (ba[offset]);
      if (block == 0 && bint == 0) {
          return null;  // Should only happen in the linear index?
      } else {
          return new Vob(block, bint);
      }
  }
  
  function unbgzf(data, lim) {
      lim = Math.min(lim || 1, data.byteLength - 50);
      var oBlockList = [];
      var ptr = [0];
      var totalSize = 0;
  
      while (ptr[0] < lim) {
          var ba = new Uint8Array(data, ptr[0], 12); // FIXME is this enough for all credible BGZF block headers?
          var xlen = (ba[11] << 8) | (ba[10]);
          // dlog('xlen[' + (ptr[0]) +']=' + xlen);
          var unc = jszlib_inflate_buffer(data, 12 + xlen + ptr[0], Math.min(65536, data.byteLength - 12 - xlen - ptr[0]), ptr);
          ptr[0] += 8;
          totalSize += unc.byteLength;
          oBlockList.push(unc);
      }
  
      if (oBlockList.length == 1) {
          return oBlockList[0];
      } else {
          var out = new Uint8Array(totalSize);
          var cursor = 0;
          for (var i = 0; i < oBlockList.length; ++i) {
              var b = new Uint8Array(oBlockList[i]);
              arrayCopy(b, 0, out, cursor, b.length);
              cursor += b.length;
          }
          return out.buffer;
      }
  }
  
  //
  // Binning (transliterated from SAM1.3 spec)
  //
  
  /* calculate bin given an alignment covering [beg,end) (zero-based, half-close-half-open) */
  function reg2bin(beg, end)
  {
      --end;
      if (beg>>14 == end>>14) return ((1<<15)-1)/7 + (beg>>14);
      if (beg>>17 == end>>17) return ((1<<12)-1)/7 + (beg>>17);
      if (beg>>20 == end>>20) return ((1<<9)-1)/7 + (beg>>20);
      if (beg>>23 == end>>23) return ((1<<6)-1)/7 + (beg>>23);
      if (beg>>26 == end>>26) return ((1<<3)-1)/7 + (beg>>26);
      return 0;
  }
  
  /* calculate the list of bins that may overlap with region [beg,end) (zero-based) */
  var MAX_BIN = (((1<<18)-1)/7);
  function reg2bins(beg, end) 
  {
      var i = 0, k, list = [];
      --end;
      list.push(0);
      for (k = 1 + (beg>>26); k <= 1 + (end>>26); ++k) list.push(k);
      for (k = 9 + (beg>>23); k <= 9 + (end>>23); ++k) list.push(k);
      for (k = 73 + (beg>>20); k <= 73 + (end>>20); ++k) list.push(k);
      for (k = 585 + (beg>>17); k <= 585 + (end>>17); ++k) list.push(k);
      for (k = 4681 + (beg>>14); k <= 4681 + (end>>14); ++k) list.push(k);
      return list;
  }