@@ -4,7 +4,7 @@ PKG_OBJECTS = \

 UCSC_OBJECTS = \

   bPlusTree.o bbiRead.o bbiWrite.o bwgCreate.o bwgQuery.o \

-  cirTree.o common.o dnaseq.o dnautil.o errabort.o hash.o linefile.o localmem.o\

+  cirTree.o common.o dnaseq.o dnautil.o errAbort.o hash.o linefile.o localmem.o\

   sqlNum.o zlibFace.o dystring.o hmmstats.o obscure.o pipeline.o \

   rangeTree.o rbTree.o memalloc.o dlist.o udc.o net.o bits.o twoBit.o \

   _cheapcgi.o internet.o https.o base64.o verbose.o os.o wildcmp.o _portimpl.o

@@ -5,7 +5,6 @@

 #include "ucsc/bbiFile.h"

 #include "ucsc/bigWig.h"

 #include "ucsc/bwgInternal.h"

-#include "ucsc/_bwgInternal.h"

 #include "bigWig.h"

 #include "handlers.h"

@@ -190,6 +189,7 @@ SEXP BWGSectionList_write(SEXP r_sections, SEXP r_seqlengths, SEXP r_compress,

   }

   pushRHandlers();

   bwgCreate(sections, lenHash, blockSize, itemsPerSlot, asLogical(r_compress),

+            FALSE /*keepAllChromosomes*/, TRUE /*fixedSummaries*/,

             (char *)CHAR(asChar(r_file)));

   freeHash(&lenHash);

   popRHandlers();

@@ -385,7 +385,7 @@ SEXP BWGFile_fromWIG(SEXP r_infile, SEXP r_seqlengths, SEXP r_outfile) {

   struct bwgSection *sections =

     bwgParseWig((char *)CHAR(asChar(r_infile)), FALSE, lenHash, itemsPerSlot,

                 lm);

-  bwgCreate(sections, lenHash, blockSize, itemsPerSlot, TRUE,

+  bwgCreate(sections, lenHash, blockSize, itemsPerSlot, TRUE, TRUE, TRUE,

             (char *)CHAR(asChar(r_outfile)));

   lmCleanup(&lm);

   freeHash(&lenHash);

@@ -1,5 +1,5 @@

 #include "ucsc/common.h"

-#include "ucsc/errabort.h"

+#include "ucsc/errAbort.h"

 #include "handlers.h"

@@ -1,3 +1,6 @@

+/* Copyright (C) 2013 The Regents of the University of California 

+ * See README in this or parent directory for licensing information. */

+

 /* bptFile - B+ Trees.  These are a method of indexing data similar to binary trees, but 

  * with many children rather than just two at each node. They work well when stored on disk,

  * since typically only two or three disk accesses are needed to locate any particular

@@ -162,7 +165,71 @@ else

     }

 }

-void rTraverse(struct bptFile *bpt, bits64 blockStart, void *context, 

+static void rFindMulti(struct bptFile *bpt, bits64 blockStart, void *key, struct slRef **pList)

+/* Find values corresponding to key and add them to pList.  You'll need to 

+ * Do a slRefFreeListAndVals() on the list when done. */

+{

+/* Seek to start of block. */

+udcSeek(bpt->udc, blockStart);

+

+/* Read block header. */

+UBYTE isLeaf;

+UBYTE reserved;

+bits16 i, childCount;

+udcMustReadOne(bpt->udc, isLeaf);

+udcMustReadOne(bpt->udc, reserved);

+boolean isSwapped = bpt->isSwapped;

+childCount = udcReadBits16(bpt->udc, isSwapped);

+

+int keySize = bpt->keySize;

+UBYTE keyBuf[keySize];   /* Place to put a key, buffered on stack. */

+UBYTE valBuf[bpt->valSize];   /* Place to put a value, buffered on stack. */

+

+if (isLeaf)

+    {

+    for (i=0; i<childCount; ++i)

+        {

+	udcMustRead(bpt->udc, keyBuf, keySize);

+	udcMustRead(bpt->udc, valBuf, bpt->valSize);

+	if (memcmp(key, keyBuf, keySize) == 0)

+	    {

+	    void *val = cloneMem(valBuf, bpt->valSize);

+	    refAdd(pList, val);

+	    }

+	}

+    }

+else

+    {

+    /* Read first key and first file offset. */

+    udcMustRead(bpt->udc, keyBuf, keySize);

+    bits64 lastFileOffset = udcReadBits64(bpt->udc, isSwapped);

+    bits64 fileOffset = lastFileOffset;

+    int lastCmp = memcmp(key, keyBuf, keySize);

+

+    /* Loop through remainder. */

+    for (i=1; i<childCount; ++i)

+	{

+	udcMustRead(bpt->udc, keyBuf, keySize);

+	fileOffset = udcReadBits64(bpt->udc, isSwapped);

+	int cmp = memcmp(key, keyBuf, keySize);

+	if (lastCmp >= 0 && cmp <= 0)

+	    {

+	    bits64 curPos = udcTell(bpt->udc);

+	    rFindMulti(bpt, lastFileOffset, key, pList);

+	    udcSeek(bpt->udc, curPos);

+	    }

+	if (cmp < 0)

+	    return;

+	lastCmp = cmp;

+	lastFileOffset = fileOffset;

+	}

+    /* If made it all the way to end, do last one too. */

+    rFindMulti(bpt, fileOffset, key, pList);

+    }

+}

+

+

+static void rTraverse(struct bptFile *bpt, bits64 blockStart, void *context, 

     void (*callback)(void *context, void *key, int keySize, void *val, int valSize) )

 /* Recursively go across tree, calling callback at leaves. */

 {

@@ -203,19 +270,83 @@ else

     }

 }

-boolean bptFileFind(struct bptFile *bpt, void *key, int keySize, void *val, int valSize)

-/* Find value associated with key.  Return TRUE if it's found. 

-*  Parameters:

-*     bpt - file handle returned by bptFileOpen

-*     key - pointer to key string, which needs to be bpt->keySize long

-*     val - pointer to where to put retrieved value

-*/

+static bits64 bptDataStart(struct bptFile *bpt)

+/* Return offset of first bit of data (as opposed to index) in file.  In hind sight I wish

+ * this were stored in the header, but fortunately it's not that hard to compute. */

+{

+bits64 offset = bpt->rootOffset;

+for (;;)

+    {

+    /* Seek to block start */

+    udcSeek(bpt->udc, offset);

+

+    /* Read block header,  break if we are leaf. */

+    UBYTE isLeaf;

+    UBYTE reserved;

+    bits16 childCount;

+    udcMustReadOne(bpt->udc, isLeaf);

+    if (isLeaf)

+         break;

+    udcMustReadOne(bpt->udc, reserved);

+    boolean isSwapped = bpt->isSwapped;

+    childCount = udcReadBits16(bpt->udc, isSwapped);

+

+    /* Read and discard first key. */

+    char keyBuf[bpt->keySize];

+    udcMustRead(bpt->udc, keyBuf, bpt->keySize);

+

+    /* Get file offset of sub-block. */

+    offset = udcReadBits64(bpt->udc, isSwapped);

+    }

+return offset;

+}

+

+static bits64 bptDataOffset(struct bptFile *bpt, bits64 itemPos)

+/* Return position of file of data corresponding to given itemPos.  For first piece of

+ * data pass in 0. */

+{

+if (itemPos >= bpt->itemCount)

+    errAbort("Item index %lld greater than item count %lld in %s", 

+	itemPos, bpt->itemCount, bpt->fileName);

+bits64 blockPos = itemPos/bpt->blockSize;

+bits32 insidePos = itemPos - blockPos*bpt->blockSize;

+int blockHeaderSize = 4;

+bits64 itemByteSize = bpt->valSize + bpt->keySize;

+bits64 blockByteSize = bpt->blockSize * itemByteSize + blockHeaderSize;

+bits64 blockOffset = blockByteSize*blockPos + bptDataStart(bpt);

+bits64 itemOffset = blockOffset + blockHeaderSize + itemByteSize * insidePos;

+return itemOffset;

+}

+

+void bptKeyAtPos(struct bptFile *bpt, bits64 itemPos, void *result)

+/* Fill in result with the key at given itemPos.  For first piece of data itemPos is 0 

+ * Result must be at least bpt->keySize.  If result is a string it won't be zero terminated

+ * by this routine.  Use bptStringKeyAtPos instead. */

+{

+bits64 offset = bptDataOffset(bpt, itemPos);

+udcSeek(bpt->udc, offset);

+udcMustRead(bpt->udc, result, bpt->keySize);

+}

+

+void bptStringKeyAtPos(struct bptFile *bpt, bits64 itemPos, char *result, int maxResultSize)

+/* Fill in result with the key at given itemPos.  The maxResultSize should be 1+bpt->keySize

+ * to accommodate zero termination of string. */

+{

+assert(maxResultSize > bpt->keySize);

+bptKeyAtPos(bpt, itemPos, result);

+result[bpt->keySize] = 0;

+}

+

+static boolean bptFileFindMaybeMulti(struct bptFile *bpt, void *key, int keySize, int valSize,

+    boolean multi, void *singleVal, struct slRef **multiVal)

+/* Do either a single or multiple find depending in multi parameter.  Only one of singleVal

+ * or multiVal should be non-NULL, depending on the same parameter. */

 {

 /* Check key size vs. file key size, and act appropriately.  If need be copy key to a local

  * buffer and zero-extend it. */

 if (keySize > bpt->keySize)

     return FALSE;

-char keyBuf[keySize];

+char keyBuf[bpt->keySize];

 if (keySize != bpt->keySize)

     {

     memcpy(keyBuf, key, keySize);

@@ -228,8 +359,34 @@ if (valSize != bpt->valSize)

     errAbort("Value size mismatch between bptFileFind (valSize=%d) and %s (valSize=%d)",

     	valSize, bpt->fileName, bpt->valSize);

-/* Call recursive finder. */

-return rFind(bpt, bpt->rootOffset, key, val);

+if (multi)

+    {

+    rFindMulti(bpt, bpt->rootOffset, key, multiVal);

+    return *multiVal != NULL;

+    }

+else

+    return rFind(bpt, bpt->rootOffset, key, singleVal);

+}

+

+boolean bptFileFind(struct bptFile *bpt, void *key, int keySize, void *val, int valSize)

+/* Find value associated with key.  Return TRUE if it's found. 

+*  Parameters:

+*     bpt - file handle returned by bptFileOpen

+*     key - pointer to key string, which needs to be bpt->keySize long

+*     val - pointer to where to put retrieved value

+*/

+{

+return bptFileFindMaybeMulti(bpt, key, keySize, valSize, FALSE, val, NULL);

+}

+

+struct slRef *bptFileFindMultiple(struct bptFile *bpt, void *key, int keySize, int valSize)

+/* Find all values associated with key.  Store this in ->val item of returned list. 

+ * Do a slRefFreeListAndVals() on list when done. */

+{

+struct slRef *list = NULL;

+bptFileFindMaybeMulti(bpt, key, keySize, valSize, TRUE, NULL, &list);

+slReverse(&list);

+return list;

 }

 void bptFileTraverse(struct bptFile *bpt, void *context,

@@ -237,7 +394,7 @@ void bptFileTraverse(struct bptFile *bpt, void *context,

 /* Traverse bPlusTree on file, calling supplied callback function at each

  * leaf item. */

 {

-  rTraverse(bpt, bpt->rootOffset, context, callback);

+return rTraverse(bpt, bpt->rootOffset, context, callback);

 }

@@ -252,10 +409,10 @@ void bptFileTraverse(struct bptFile *bpt, void *context,

  *  01 02 03 04   05 06 07 08  09 10 11 12   13 14 15 16   17 18 19 20   21 22 23 24  25 26 27

  */

-static int xToY(int x, unsigned y)

+static long xToY(int x, unsigned y)

 /* Return x to the Y power, with y usually small. */

 {

-int i, val = 1;

+long i, val = 1;

 for (i=0; i<y; ++i)

     val *= x;

 return val;

@@ -275,7 +432,7 @@ return levels;

 static bits32 writeIndexLevel(bits16 blockSize, 

-	void *itemArray, int itemSize, int itemCount, 

+	void *itemArray, int itemSize, long itemCount, 

 	bits32 indexOffset, int level, 

 	void (*fetchKey)(const void *va, char *keyBuf), bits32 keySize, bits32 valSize,

 	FILE *f)

@@ -284,39 +441,42 @@ static bits32 writeIndexLevel(bits16 blockSize,

 char *items = itemArray;

 /* Calculate number of nodes to write at this level. */

-int slotSizePer = xToY(blockSize, level);   // Number of items per slot in node

-int nodeSizePer = slotSizePer * blockSize;  // Number of items per node

-int nodeCount = (itemCount + nodeSizePer - 1)/nodeSizePer;	

+long slotSizePer = xToY(blockSize, level);   // Number of items per slot in node

+long nodeSizePer = slotSizePer * blockSize;  // Number of items per node

+long nodeCount = (itemCount + nodeSizePer - 1)/nodeSizePer;	

+

 /* Calculate sizes and offsets. */

-int bytesInIndexBlock = (bptBlockHeaderSize + blockSize * (keySize+sizeof(bits64)));

-int bytesInLeafBlock = (bptBlockHeaderSize + blockSize * (keySize+valSize));

+long bytesInIndexBlock = (bptBlockHeaderSize + blockSize * (keySize+sizeof(bits64)));

+long bytesInLeafBlock = (bptBlockHeaderSize + blockSize * (keySize+valSize));

 bits64 bytesInNextLevelBlock = (level == 1 ? bytesInLeafBlock : bytesInIndexBlock);

 bits64 levelSize = nodeCount * bytesInIndexBlock;

 bits64 endLevel = indexOffset + levelSize;

 bits64 nextChild = endLevel;

+

 UBYTE isLeaf = FALSE;

 UBYTE reserved = 0;

-int i,j;

+long i,j;

 char keyBuf[keySize+1];

 keyBuf[keySize] = 0;

 for (i=0; i<itemCount; i += nodeSizePer)

     {

     /* Calculate size of this block */

-    bits16 countOne = (itemCount - i + slotSizePer - 1)/slotSizePer;

+    long countOne = (itemCount - i + slotSizePer - 1)/slotSizePer;

     if (countOne > blockSize)

         countOne = blockSize;

+    bits16 shortCountOne = countOne;

     /* Write block header. */

     writeOne(f, isLeaf);

     writeOne(f, reserved);

-    writeOne(f, countOne);

+    writeOne(f, shortCountOne);

     /* Write out the slots that are used one by one, and do sanity check. */

     int slotsUsed = 0;

-    int endIx = i + nodeSizePer;

+    long endIx = i + nodeSizePer;

     if (endIx > itemCount)

         endIx = itemCount;

     for (j=i; j<endIx; j += slotSizePer)

@@ -329,7 +489,7 @@ for (i=0; i<itemCount; i += nodeSizePer)

 	nextChild += bytesInNextLevelBlock;

 	++slotsUsed;

 	}

-    assert(slotsUsed == countOne);

+    assert(slotsUsed == shortCountOne);

     /* Write out empty slots as all zero. */

     int slotSize = keySize + sizeof(bits64);

@@ -76,19 +76,32 @@ void bptFileDetach(struct bptFile **pBpt);

 boolean bptFileFind(struct bptFile *bpt, void *key, int keySize, void *val, int valSize);

 /* Find value associated with key.  Return TRUE if it's found. 

-*  Parameters:

-*     bpt - file handle returned by bptFileOpen

-*     key - pointer to key string

-*     keySize - size of key.  Normally just strlen(key)

-*     val - pointer to where to put retrieved value

-*     valSize - size of memory buffer that will hold val.  Should match bpt->valSize.

-*/

+ *  Parameters:

+ *     bpt - file handle returned by bptFileOpen

+ *     key - pointer to key string

+ *     keySize - size of key.  Normally just strlen(key)

+ *     val - pointer to where to put retrieved value

+ *     valSize - size of memory buffer that will hold val.  Should match bpt->valSize.

+ */

+

+struct slRef *bptFileFindMultiple(struct bptFile *bpt, void *key, int keySize, int valSize);

+/* Find all values associated with key.  Store this in ->val item of returned list. 

+ * Do a slRefFreeListAndVals() on list when done. */

 void bptFileTraverse(struct bptFile *bpt, void *context,

     void (*callback)(void *context, void *key, int keySize, void *val, int valSize) );

 /* Traverse bPlusTree on file, calling supplied callback function at each

  * leaf item. */

+void bptKeyAtPos(struct bptFile *bpt, bits64 itemPos, void *result);

+/* Fill in result with the key at given itemPos.  For first piece of data itemPos is 0 

+ * and for last piece is bpt->itemCount - 1.  Result must be at least bpt->keySize.  

+ * If result is a string it won't be zero terminated

+ * by this routine.  Use bptStringKeyAtPos instead. */

+

+void bptStringKeyAtPos(struct bptFile *bpt, bits64 itemPos, char *result, int maxResultSize);

+/* Fill in result with the key at given itemPos.  The maxResultSize should be 1+bpt->keySize

+ * to accommodate zero termination of string. */

 void bptFileCreate(

 	void *itemArray, 	/* Sorted array of things to index. */

@@ -1,7 +1,9 @@

+/* Copyright (C) 2011 The Regents of the University of California 

+ * See README in this or parent directory for licensing information. */

+

 #include "common.h"

 #include "base64.h"

-static char const rcsid[] = "$Id: base64.c,v 1.6 2008/09/17 18:00:47 galt Exp $";

 char *base64Encode(char *input, size_t inplen)

 /* Use base64 to encode a string.  Returns one long encoded

@@ -4,6 +4,8 @@

 #define BBIFILE_H

 #include "cirTree.h"

+#include "linefile.h"

+#include "localmem.h"

 /* bigWig/bigBed file structure:

  *     fixedWidthHeader

@@ -18,7 +20,7 @@

  *         autoSqlOffset        8 bytes (for bigWig 0) (0 if no autoSql information)

  *         totalSummaryOffset   8 bytes (0 in earlier versions of file lacking totalSummary)

  *         uncompressBufSize    4 bytes (Size of uncompression buffer.  0 if uncompressed.)

- *         reserved             8 bytes (0 for now)

+ *         extensionOffset      8 bytes (Offset to header extension 0 if no such extension)

  *     zoomHeaders		there are zoomLevels number of these

  *         reductionLevel	4 bytes

  *	   reserved		4 bytes

@@ -31,6 +33,19 @@

  *         maxVal              8 bytes float (for bigBed maximum depth of coverage)

  *         sumData             8 bytes float (for bigBed sum of coverage)

  *         sumSquared          8 bytes float (for bigBed sum of coverage squared)

+ *     extendedHeader

+ *         extensionSize       2 size of extended header in bytes - currently 64

+ *         extraIndexCount     2 number of extra fields we will be indexing

+ *         extraIndexListOffset 8 Offset to list of non-chrom/start/end indexes

+ *         reserved            48 All zeroes for now

+ *     extraIndexList - one of these for each extraIndex 

+ *         type                2 Type of index.  Always 0 for bPlusTree now

+ *         fieldCount          2 Number of fields used in this index.  Always 1 for now

+ *         indexOffset         8 offset for this index in file

+ *         reserved            4 All zeroes for now

+ *         fieldList - one of these for each field being used in _this_ index

+ *            fieldId          2 index of field within record

+ *            reserved         2 All zeroes for now

  *     chromosome b+ tree       bPlusTree index

  *     full data

  *         sectionCount		8 bytes (item count for bigBeds)

@@ -49,6 +64,7 @@

  *                 sumData      4 bytes float

  *                 sumSquares   4 bytes float

  *         zoom index        	cirTree index

+ *     extraIndexes [optional]  bPlusTreeIndex for each extra field that is indexed

  *     magic# 		4 bytes - same as magic number at start of header

  */

@@ -102,10 +118,17 @@ struct bbiFile

     bits16 fieldCount;		/* Number of columns in bed version. */

     bits16 definedFieldCount;   /* Number of columns using bed standard definitions. */

     bits64 asOffset;		/* Offset to embedded null-terminated AutoSQL file. */

-    bits64 totalSummaryOffset;	/* Offset to total summary information if any.  (On older files have to calculate) */

+    bits64 totalSummaryOffset;	/* Offset to total summary information if any.  

+				   (On older files have to calculate) */

     bits32 uncompressBufSize;	/* Size of uncompression buffer, 0 if uncompressed */

+    bits64 extensionOffset;	/* Start of header extension block or 0 if none. */

     struct cirTreeFile *unzoomedCir;	/* Unzoomed data index in memory - may be NULL. */

     struct bbiZoomLevel *levelList;	/* List of zoom levels. */

+

+    /* Fields based on extension block. */

+    bits16 extensionSize;   /* Size of extension block */

+    bits16 extraIndexCount; /* Number of extra indexes (on fields other than chrom,start,end */ 

+    bits64 extraIndexListOffset;    /* Offset to list of extra indexes */

     };

@@ -150,6 +173,14 @@ void bbiChromInfoKey(const void *va, char *keyBuf);

 void *bbiChromInfoVal(const void *va);

 /* Get val field out of bbiChromInfo. */

+char *bbiCachedChromLookup(struct bbiFile *bbi, int chromId, int lastChromId,

+    char *chromBuf, int chromBufSize);

+/* Return chromosome name corresponding to chromId.  Because this is a bit expensive,

+ * if you are doing this repeatedly pass in the chromId used in the previous call to

+ * this in lastChromId,  which will save it from doing the lookup again on the same

+ * chromosome.  Pass in -1 to lastChromId if this is the first time or if you can't be

+ * bothered.  The chromBufSize should be at greater or equal to bbi->keySize+1.  */

+

 struct bbiChromUsage

 /* Information on how many items per chromosome etc.  Used by multipass bbiFile writers. */

     {

@@ -249,7 +280,7 @@ boolean bbiSummaryArray(struct bbiFile *bbi, char *chrom, bits32 start, bits32 e

 struct bbiSummaryElement bbiTotalSummary(struct bbiFile *bbi);

 /* Return summary of entire file! */

-/****** Write side of things - implemented in bbiWrite.c ********/

+/****** Write side of things - implemented in bbiWrite.c.  Few people need this. ********/

 struct bbiBoundsArray

 /* Minimum info needed for r-tree indexer - where a section lives on disk and the

@@ -287,7 +318,7 @@ void bbiSumOutStreamWrite(struct bbiSumOutStream *stream, struct bbiSummary *sum

 void bbiOutputOneSummaryFurtherReduce(struct bbiSummary *sum, 

 	struct bbiSummary **pTwiceReducedList, 

 	int doubleReductionSize, struct bbiBoundsArray **pBoundsPt, 

-	struct bbiBoundsArray *boundsEnd, bits32 chromSize, struct lm *lm, 

+	struct bbiBoundsArray *boundsEnd, struct lm *lm, 

 	struct bbiSumOutStream *stream);

 /* Write out sum to file, keeping track of minimal info on it in *pBoundsPt, and also adding

  * it to second level summary. */

@@ -296,8 +327,6 @@ struct bbiSummary *bbiSummarySimpleReduce(struct bbiSummary *list, int reduction

 /* Do a simple reduction - where among other things the reduction level is an integral

  * multiple of the previous reduction level, and the list is sorted. Allocate result out of lm. */

-#define bbiMaxZoomLevels 10	/* Maximum zoom levels produced by writers. */

-

 void bbiWriteDummyHeader(FILE *f);

 /* Write out all-zero header, just to reserve space for it. */

@@ -325,9 +354,69 @@ void bbiChromUsageFree(struct bbiChromUsage **pUsage);

 void bbiChromUsageFreeList(struct bbiChromUsage **pList);

 /* free a list of bbiChromUsage structures */

-struct bbiChromUsage *bbiChromUsageFromBedFile(struct lineFile *lf, 

-	struct hash *chromSizesHash, int *retMinDiff, double *retAveSize, bits64 *retBedCount);

-/* Go through bed file and collect chromosomes and statistics. Free with bbiChromUsageFreeList */

+struct bbNamedFileChunk 

+/* A name associated with an offset into a possibly large file.  Used for extra

+ * indexes in bigBed files. */

+    {

+    char *name;	    /* Name of chunk. */

+    bits64 offset;  /* Start in file. */

+    bits64 size;    /* Size in file. */

+    };

+

+struct bbExIndexMaker

+/* A helper structure to make indexes beyond primary one.  Just used for bigBeds */

+    {

+    bits16 indexCount;          /* Number of extra indexes. */

+        /* Kind of wish next four fields,  all of which are arrays indexed

+         * by the same thing,  were a single array of a structure instead. */

+    bits16 *indexFields;        /* array of field ids, one for each extra index. */

+    int *maxFieldSize;          /* array of maximum sizes seen for this field. */

+    struct bbNamedFileChunk **chunkArrayArray; /* where we keep name/start/size triples */

+    bits64 *fileOffsets;        /* array of file offsets where indexes starts. */

+    int recordCount;            /* number of records in file. */

+    };

+

+struct bbiChromUsage *bbiChromUsageFromBedFile(struct lineFile *lf, struct hash *chromSizesHash, 

+	struct bbExIndexMaker *eim, int *retMinDiff, double *retAveSize, bits64 *retBedCount);

+/* Go through bed file and collect chromosomes and statistics.  If eim parameter is non-NULL

+ * collect max field sizes there too. */

+

+#define bbiMaxZoomLevels 10	/* Max number of zoom levels */

+#define bbiResIncrement 4	/* Amount to reduce at each zoom level */

+

+int bbiCalcResScalesAndSizes(int aveSize, 

+    int resScales[bbiMaxZoomLevels], int resSizes[bbiMaxZoomLevels]);

+/* Fill in resScales with amount to zoom at each level, and zero out resSizes based

+ * on average span. Returns the number of zoom levels we actually will use. */

+

+typedef struct bbiSummary *bbiWriteReducedOnceReturnReducedTwice(

+	struct bbiChromUsage *usageList, int fieldCount,

+	struct lineFile *lf, bits32 initialReduction, bits32 initialReductionCount, 

+	int zoomIncrement, int blockSize, int itemsPerSlot, boolean doCompress,

+	struct lm *lm, FILE *f, bits64 *retDataStart, bits64 *retIndexStart,

+	struct bbiSummaryElement *totalSum);

+/* Typedef for a function that writes out data reduced by factor of initial reduction, and

+ * also returns an array of bbiSummaries for the next reduction level. */

+

+int bbiWriteZoomLevels(

+    struct lineFile *lf,    /* Input file. */

+    FILE *f,		    /* Output. */

+    int blockSize,	    /* Size of index block */

+    int itemsPerSlot,	    /* Number of data points bundled at lowest level. */

+    bbiWriteReducedOnceReturnReducedTwice writeReducedOnceReturnReducedTwice,   /* callback */

+    int fieldCount,	    /* Number of fields in bed (4 for bedGraph) */

+    boolean doCompress,	    /* Do we compress.  Answer really should be yes! */

+    bits64 dataSize,	    /* Size of data on disk (after compression if any). */

+    struct bbiChromUsage *usageList, /* Result from bbiChromUsageFromBedFile */

+    int resTryCount, int resScales[], int resSizes[],   /* How much to zoom at each level */

+    bits32 zoomAmounts[bbiMaxZoomLevels],      /* Fills in amount zoomed at each level. */

+    bits64 zoomDataOffsets[bbiMaxZoomLevels],  /* Fills in where data starts for each zoom level. */

+    bits64 zoomIndexOffsets[bbiMaxZoomLevels], /* Fills in where index starts for each level. */

+    struct bbiSummaryElement *totalSum);

+/* Write out all the zoom levels and return the number of levels written.  Writes 

+ * actual zoom amount and the offsets of the zoomed data and index in the last three

+ * parameters.  Sorry for all the parameters - it was this or duplicate a big chunk of

+ * code between bedToBigBed and bedGraphToBigWig. */

 int bbiCountSectionsNeeded(struct bbiChromUsage *usageList, int itemsPerSlot);

 /* Count up number of sections needed for data. */

@@ -358,4 +447,8 @@ boolean bbiFileCheckSigs(char *fileName, bits32 sig, char *typeName);

 time_t bbiUpdateTime(struct bbiFile *bbi);

 /* return bbi->udc->updateTime */

+struct bbiSummary *bbiSummariesInRegion(struct bbiZoomLevel *zoom, struct bbiFile *bbi,

+        int chromId, bits32 start, bits32 end);

+/* Return list of all summaries in region at given zoom level of bbiFile. */

+

 #endif /* BBIFILE_H */

@@ -1,11 +1,13 @@

 /* bbiRead - Big Binary Indexed file.  Stuff that's common between bigWig and bigBed on the 

  * read side. */

+/* Copyright (C) 2014 The Regents of the University of California 

+ * See README in this or parent directory for licensing information. */

+

 #include "common.h"

 #include "linefile.h"

 #include "hash.h"

 #include "obscure.h"

-#include "localmem.h"

 #include "zlibFace.h"

 #include "bPlusTree.h"

 #include "hmmstats.h"

@@ -115,9 +117,7 @@ bbi->definedFieldCount = udcReadBits16(udc, isSwapped);

 bbi->asOffset = udcReadBits64(udc, isSwapped);

 bbi->totalSummaryOffset = udcReadBits64(udc, isSwapped);

 bbi->uncompressBufSize = udcReadBits32(udc, isSwapped);

-

-/* Skip over reserved area. */

-udcSeek(udc, 64);

+bbi->extensionOffset = udcReadBits64(udc, isSwapped);

 /* Read zoom headers. */

 int i;

@@ -134,6 +134,15 @@ for (i=0; i<bbi->zoomLevels; ++i)

 slReverse(&levelList);

 bbi->levelList = levelList;

+/* Deal with header extension if any. */

+if (bbi->extensionOffset != 0)

+    {

+    udcSeek(udc, bbi->extensionOffset);

+    bbi->extensionSize = udcReadBits16(udc, isSwapped);

+    bbi->extraIndexCount = udcReadBits16(udc, isSwapped);

+    bbi->extraIndexListOffset = udcReadBits64(udc, isSwapped);

+    }

+

 /* Attach B+ tree of chromosome names and ids. */

 udcSeek(udc, bbi->chromTreeOffset);

 bbi->chromBpt =  bptFileAttach(fileName, udc);

@@ -157,6 +166,16 @@ if (bwf != NULL)

     }

 }

+static void chromIdSizeHandleSwapped(boolean isSwapped, struct bbiChromIdSize *idSize)

+/* Swap bytes in chromosome Id and Size as needed. */

+{

+if (isSwapped)

+    {

+    idSize->chromId = byteSwap32(idSize->chromId);

+    idSize->chromSize = byteSwap32(idSize->chromSize);

+    }

+}

+

 struct fileOffsetSize *bbiOverlappingBlocks(struct bbiFile *bbi, struct cirTreeFile *ctf,

 	char *chrom, bits32 start, bits32 end, bits32 *retChromId)

@@ -165,8 +184,7 @@ struct fileOffsetSize *bbiOverlappingBlocks(struct bbiFile *bbi, struct cirTreeF

 struct bbiChromIdSize idSize;

 if (!bptFileFind(bbi->chromBpt, chrom, strlen(chrom), &idSize, sizeof(idSize)))

     return NULL;

-if (bbi->isSwapped)

-    idSize.chromId = byteSwap32(idSize.chromId);

+chromIdSizeHandleSwapped(bbi->isSwapped, &idSize);

 if (retChromId != NULL)

     *retChromId = idSize.chromId;

 return cirTreeFindOverlappingBlocks(ctf, idSize.chromId, start, end);

@@ -186,15 +204,11 @@ struct chromNameCallbackContext *c = context;

 struct bbiChromInfo *info;

 struct bbiChromIdSize *idSize = val;

 assert(valSize == sizeof(*idSize));

+chromIdSizeHandleSwapped(c->isSwapped, idSize);

 AllocVar(info);

 info->name = cloneStringZ(key, keySize);

 info->id = idSize->chromId;

 info->size = idSize->chromSize;

-if (c->isSwapped)

-    {

-    info->id = byteSwap32(info->id);

-    info->size = byteSwap32(info->size);

-    }

 slAddHead(&c->list, info);

 }

@@ -215,6 +229,7 @@ bits32 bbiChromSize(struct bbiFile *bbi, char *chrom)

 struct bbiChromIdSize idSize;

 if (!bptFileFind(bbi->chromBpt, chrom, strlen(chrom), &idSize, sizeof(idSize)))

     return 0;

+chromIdSizeHandleSwapped(bbi->isSwapped, &idSize);

 return idSize.chromSize;

 }

@@ -303,6 +318,7 @@ sum->chromId = udcReadBits32(udc, isSwapped);

 sum->start = udcReadBits32(udc, isSwapped);

 sum->end = udcReadBits32(udc, isSwapped);

 sum->validCount = udcReadBits32(udc, isSwapped);

+// looks like a bug to me, these should call udcReadFloat() 

 udcMustReadOne(udc, sum->minVal);

 udcMustReadOne(udc, sum->maxVal);

 udcMustReadOne(udc, sum->sumData);

@@ -310,6 +326,22 @@ udcMustReadOne(udc, sum->sumSquares);

 }

 #endif /* UNUSED */

+static void bbiSummaryHandleSwapped(struct bbiFile *bbi, struct bbiSummaryOnDisk *in)

+/* Swap integer fields in summary as needed. */

+{

+if (bbi->isSwapped)

+    {

+    in->chromId = byteSwap32(in->chromId);

+    in->start = byteSwap32(in->start);

+    in->end = byteSwap32(in->end);

+    in->validCount = byteSwap32(in->validCount);

+    in->minVal = byteSwapFloat(in->minVal);

+    in->maxVal = byteSwapFloat(in->maxVal);

+    in->sumData = byteSwapFloat(in->sumData);

+    in->sumSquares = byteSwapFloat(in->sumSquares);

+    }

+}

+

 static struct bbiSummary *bbiSummaryFromOnDisk(struct bbiSummaryOnDisk *in)

 /* Create a bbiSummary unlinked to anything from input in onDisk format. */

 {

@@ -326,7 +358,7 @@ out->sumSquares = in->sumSquares;

 return out;

 }

-static struct bbiSummary *bbiSummariesInRegion(struct bbiZoomLevel *zoom, struct bbiFile *bbi, 

+struct bbiSummary *bbiSummariesInRegion(struct bbiZoomLevel *zoom, struct bbiFile *bbi, 

 	int chromId, bits32 start, bits32 end)

 /* Return list of all summaries in region at given zoom level of bbiFile. */

 {

@@ -386,6 +418,7 @@ for (block = blockList; block != NULL; )

 	    {

 	    dSum = (void *)blockPt;

 	    blockPt += sizeof(*dSum);

+	    bbiSummaryHandleSwapped(bbi, dSum);

 	    if (dSum->chromId == chromId)

 		{

 		int s = max(dSum->start, start);

@@ -409,6 +442,23 @@ slReverse(&sumList);

 return sumList;

 }

+static int normalizeCount(struct bbiSummaryElement *el, double countFactor, 

+    double minVal, double maxVal, double sumData, double sumSquares)

+/* normalize statistics to be based on an integer number of valid bases 

+ * Integer value is the smallest integer not less than countFactor */

+{

+bits32 validCount = ceil(countFactor);

+double normFactor = (double)validCount/countFactor;

+

+el->validCount = validCount;

+el->minVal = minVal;

+el->maxVal = maxVal;

+el->sumData = sumData * normFactor;

+el->sumSquares = sumSquares * normFactor;

+

+return validCount;

+}

+

 static bits32 bbiSummarySlice(struct bbiFile *bbi, bits32 baseStart, bits32 baseEnd, 

 	struct bbiSummary *sumList, struct bbiSummaryElement *el)

 /* Update retVal with the average value if there is any data in interval.  Return number

@@ -421,6 +471,7 @@ if (sumList != NULL)

     double minVal = sumList->minVal;

     double maxVal = sumList->maxVal;

     double sumData = 0, sumSquares = 0;

+    double countFactor = 0.0;

     struct bbiSummary *sum;

     for (sum = sumList; sum != NULL && sum->start < baseEnd; sum = sum->next)

@@ -429,7 +480,7 @@ if (sumList != NULL)

 	if (overlap > 0)

 	    {

 	    double overlapFactor = (double)overlap / (sum->end - sum->start);

-	    validCount += sum->validCount * overlapFactor;

+	    countFactor += sum->validCount * overlapFactor;

 	    sumData += sum->sumData * overlapFactor;

 	    sumSquares += sum->sumSquares * overlapFactor;

 	    if (maxVal < sum->maxVal)

@@ -438,24 +489,20 @@ if (sumList != NULL)

 		minVal = sum->minVal;

 	    }

 	}

-    if (validCount > 0)

-	{

-	el->validCount = validCount;

-	el->minVal = minVal;

-	el->maxVal = maxVal;

-	el->sumData = sumData;

-	el->sumSquares = sumSquares;

-	}

+

+    if (countFactor > 0)

+	validCount = normalizeCount(el, countFactor, minVal, maxVal, sumData, sumSquares);

     }

 return validCount;

 }

 static int bbiChromId(struct bbiFile *bbi, char *chrom)

-/* Return chromosome size */

+/* Return chromosome Id */

 {

 struct bbiChromIdSize idSize;

 if (!bptFileFind(bbi->chromBpt, chrom, strlen(chrom), &idSize, sizeof(idSize)))

     return -1;

+chromIdSizeHandleSwapped(bbi->isSwapped, &idSize);

 return idSize.chromId;

 }

@@ -501,10 +548,11 @@ static bits32 bbiIntervalSlice(struct bbiFile *bbi, bits32 baseStart, bits32 bas

 /* Update retVal with the average value if there is any data in interval.  Return number

  * of valid data bases in interval. */

 {

-double validCount = 0;

+bits32 validCount = 0;

 if (intervalList != NULL)

     {

+    double countFactor = 0;

     struct bbiInterval *interval;

     double sumData = 0, sumSquares = 0;

     double minVal = intervalList->val;

@@ -519,7 +567,7 @@ if (intervalList != NULL)

 	    int intervalSize = interval->end - interval->start;

 	    double overlapFactor = (double)overlap / intervalSize;

 	    double intervalWeight = intervalSize * overlapFactor;

-	    validCount += intervalWeight;

+	    countFactor += intervalWeight;

 	    sumData += interval->val * intervalWeight;

 	    sumSquares += interval->val * interval->val * intervalWeight;

 	    if (maxVal < interval->val)

@@ -528,13 +576,10 @@ if (intervalList != NULL)

 		minVal = interval->val;

 	    }

 	}

-    el->validCount = round(validCount);

-    el->minVal = minVal;

-    el->maxVal = maxVal;

-    el->sumData = sumData;

-    el->sumSquares = sumSquares;

+

+    validCount = normalizeCount(el, countFactor, minVal, maxVal, sumData, sumSquares);

     }

-return round(validCount);

+return validCount;

 }

@@ -547,7 +592,7 @@ struct bbiInterval *intervalList = NULL, *interval;

 struct lm *lm = lmInit(0);

 intervalList = (*fetchIntervals)(bbi, chrom, start, end, lm);

 boolean result = FALSE;

-if (intervalList != NULL);

+if (intervalList != NULL)

     {

     int i;

     bits32 baseStart = start, baseEnd;

@@ -736,3 +781,17 @@ time_t bbiUpdateTime(struct bbiFile *bbi)

 struct udcFile *udc = bbi->udc;

 return udcUpdateTime(udc);

 }

+

+char *bbiCachedChromLookup(struct bbiFile *bbi, int chromId, int lastChromId,

+    char *chromBuf, int chromBufSize)

+/* Return chromosome name corresponding to chromId.  Because this is a bit expensive,

+ * if you are doing this repeatedly pass in the chromId used in the previous call to

+ * this in lastChromId,  which will save it from doing the lookup again on the same

+ * chromosome.  Pass in -1 to lastChromId if this is the first time or if you can't be

+ * bothered.  The chromBufSize should be at greater or equal to bbi->keySize+1.  */

+{

+if (chromId != lastChromId)

+    bptStringKeyAtPos(bbi->chromBpt, chromId, chromBuf, chromBufSize);

+return chromBuf;

+}

+

@@ -1,8 +1,12 @@

+/* bbiWrite.c - Routines to help write bigWig and bigBed files. See also bbiFile.h */

+

+/* Copyright (C) 2014 The Regents of the University of California 

+ * See README in this or parent directory for licensing information. */