c_cpp/lib/htslib/cram/cram_stats.c

   1 /*
   2 Copyright (c) 2012-2013 Genome Research Ltd.
   3 Author: James Bonfield <jkb@sanger.ac.uk>
   4
   5 Redistribution and use in source and binary forms, with or without
   6 modification, are permitted provided that the following conditions are met:
   7
   8    1. Redistributions of source code must retain the above copyright notice,
   9 this list of conditions and the following disclaimer.
  10
  11    2. Redistributions in binary form must reproduce the above copyright notice,
  12 this list of conditions and the following disclaimer in the documentation
  13 and/or other materials provided with the distribution.
  14
  15    3. Neither the names Genome Research Ltd and Wellcome Trust Sanger
  16 Institute nor the names of its contributors may be used to endorse or promote
  17 products derived from this software without specific prior written permission.
  18
  19 THIS SOFTWARE IS PROVIDED BY GENOME RESEARCH LTD AND CONTRIBUTORS "AS IS" AND
  20 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  21 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  22 DISCLAIMED. IN NO EVENT SHALL GENOME RESEARCH LTD OR CONTRIBUTORS BE LIABLE
  23 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  24 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  25 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  26 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  27 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  28 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  29 */
  30
  31 #include <config.h>
  32
  33 #include <stdio.h>
  34 #include <errno.h>
  35 #include <assert.h>
  36 #include <stdlib.h>
  37 #include <string.h>
  38 #include <zlib.h>
  39 #include <sys/types.h>
  40 #include <sys/stat.h>
  41 #include <math.h>
  42
  43 #include "cram/cram.h"
  44 #include "cram/os.h"
  45
  46 cram_stats *cram_stats_create(void) {
  47     return calloc(1, sizeof(cram_stats));
  48 }
  49
  50 void cram_stats_add(cram_stats *st, int32_t val) {
  51     st->nsamp++;
  52
  53     //assert(val >= 0);
  54
  55     if (val < MAX_STAT_VAL && val >= 0) {
  56         st->freqs[val]++;
  57     } else {
  58         khint_t k;
  59         int r;
  60
  61         if (!st->h) {
  62             st->h = kh_init(m_i2i);
  63         }
  64
  65         k = kh_put(m_i2i, st->h, val, &r);
  66         if (r == 0)
  67             kh_val(st->h, k)++;
  68         else if (r != -1)
  69             kh_val(st->h, k) = 1;
  70         else
  71             ; // FIXME: handle error
  72     }
  73 }
  74
  75 void cram_stats_del(cram_stats *st, int32_t val) {
  76     st->nsamp--;
  77
  78     //assert(val >= 0);
  79
  80     if (val < MAX_STAT_VAL && val >= 0) {
  81         st->freqs[val]--;
  82         assert(st->freqs[val] >= 0);
  83     } else if (st->h) {
  84         khint_t k = kh_get(m_i2i, st->h, val);
  85
  86         if (k != kh_end(st->h)) {
  87             if (--kh_val(st->h, k) == 0)
  88                 kh_del(m_i2i, st->h, k);
  89         } else {
  90             fprintf(stderr, "Failed to remove val %d from cram_stats\n", val);
  91             st->nsamp++;
  92         }
  93     } else {
  94         fprintf(stderr, "Failed to remove val %d from cram_stats\n", val);
  95         st->nsamp++;
  96     }
  97 }
  98
  99 void cram_stats_dump(cram_stats *st) {
 100     int i;
 101     fprintf(stderr, "cram_stats:\n");
 102     for (i = 0; i < MAX_STAT_VAL; i++) {
 103         if (!st->freqs[i])
 104             continue;
 105         fprintf(stderr, "\t%d\t%d\n", i, st->freqs[i]);
 106     }
 107     if (st->h) {
 108         khint_t k;
 109         for (k = kh_begin(st->h); k != kh_end(st->h); k++) {
 110             if (!kh_exist(st->h, k))
 111                 continue;
 112
 113             fprintf(stderr, "\t%d\t%d\n", kh_key(st->h, k), kh_val(st->h, k));
 114         }
 115     }
 116 }
 117
 118 /*
 119  * Computes entropy from integer frequencies for various encoding methods and
 120  * picks the best encoding.
 121  *
 122  * FIXME: we could reuse some of the code here for the actual encoding
 123  * parameters too. Eg the best 'k' for SUBEXP or the code lengths for huffman.
 124  *
 125  * Returns the best codec to use.
 126  */
 127 enum cram_encoding cram_stats_encoding(cram_fd *fd, cram_stats *st) {
 128     int nvals, i, ntot = 0, max_val = 0, min_val = INT_MAX;
 129     int *vals = NULL, *freqs = NULL, vals_alloc = 0;
 130
 131     //cram_stats_dump(st);
 132
 133     /* Count number of unique symbols */
 134     for (nvals = i = 0; i < MAX_STAT_VAL; i++) {
 135         if (!st->freqs[i])
 136             continue;
 137         if (nvals >= vals_alloc) {
 138             vals_alloc = vals_alloc ? vals_alloc*2 : 1024;
 139             vals  = realloc(vals,  vals_alloc * sizeof(int));
 140             freqs = realloc(freqs, vals_alloc * sizeof(int));
 141             if (!vals || !freqs) {
 142                 if (vals)  free(vals);
 143                 if (freqs) free(freqs);
 144                 return E_HUFFMAN; // Cannot do much else atm
 145             }
 146         }
 147         vals[nvals] = i;
 148         freqs[nvals] = st->freqs[i];
 149         ntot += freqs[nvals];
 150         if (max_val < i) max_val = i;
 151         if (min_val > i) min_val = i;
 152         nvals++;
 153     }
 154     if (st->h) {
 155         khint_t k;
 156         int i;
 157
 158         for (k = kh_begin(st->h); k != kh_end(st->h); k++) {
 159             if (!kh_exist(st->h, k))
 160                 continue;
 161
 162             if (nvals >= vals_alloc) {
 163                 vals_alloc = vals_alloc ? vals_alloc*2 : 1024;
 164                 vals  = realloc(vals,  vals_alloc * sizeof(int));
 165                 freqs = realloc(freqs, vals_alloc * sizeof(int));
 166                 if (!vals || !freqs)
 167                     return E_HUFFMAN; // Cannot do much else atm
 168             }
 169             i = kh_key(st->h, k);
 170             vals[nvals]=i;
 171             freqs[nvals] = kh_val(st->h, k);
 172             ntot += freqs[nvals];
 173             if (max_val < i) max_val = i;
 174             if (min_val > i) min_val = i;
 175             nvals++;
 176         }
 177     }
 178
 179     st->nvals = nvals;
 180     assert(ntot == st->nsamp);
 181
 182     free(vals);
 183     free(freqs);
 184
 185     /*
 186      * Simple policy that everything is external unless it can be
 187      * encoded using zero bits as a unary item huffman table.
 188      */
 189     return nvals <= 1 ? E_HUFFMAN : E_EXTERNAL;
 190 }
 191
 192 void cram_stats_free(cram_stats *st) {
 193     if (st->h)
 194         kh_destroy(m_i2i, st->h);
 195     free(st);
 196 }