New upstream release.
[samtools.git] / bam_lpileup.c
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <assert.h>
4 #include "bam.h"
5 #include "ksort.h"
6
7 #define TV_GAP 2
8
9 typedef struct __freenode_t {
10         uint32_t level:28, cnt:4;
11         struct __freenode_t *next;
12 } freenode_t, *freenode_p;
13
14 #define freenode_lt(a,b) ((a)->cnt < (b)->cnt || ((a)->cnt == (b)->cnt && (a)->level < (b)->level))
15 KSORT_INIT(node, freenode_p, freenode_lt)
16
17 /* Memory pool, similar to the one in bam_pileup.c */
18 typedef struct {
19         int cnt, n, max;
20         freenode_t **buf;
21 } mempool_t;
22
23 static mempool_t *mp_init()
24 {
25         return (mempool_t*)calloc(1, sizeof(mempool_t));
26 }
27 static void mp_destroy(mempool_t *mp)
28 {
29         int k;
30         for (k = 0; k < mp->n; ++k) free(mp->buf[k]);
31         free(mp->buf); free(mp);
32 }
33 static inline freenode_t *mp_alloc(mempool_t *mp)
34 {
35         ++mp->cnt;
36         if (mp->n == 0) return (freenode_t*)calloc(1, sizeof(freenode_t));
37         else return mp->buf[--mp->n];
38 }
39 static inline void mp_free(mempool_t *mp, freenode_t *p)
40 {
41         --mp->cnt; p->next = 0; p->cnt = TV_GAP;
42         if (mp->n == mp->max) {
43                 mp->max = mp->max? mp->max<<1 : 256;
44                 mp->buf = (freenode_t**)realloc(mp->buf, sizeof(freenode_t*) * mp->max);
45         }
46         mp->buf[mp->n++] = p;
47 }
48
49 /* core part */
50 struct __bam_lplbuf_t {
51         int max, n_cur, n_pre;
52         int max_level, *cur_level, *pre_level;
53         mempool_t *mp;
54         freenode_t **aux, *head, *tail;
55         int n_nodes, m_aux;
56         bam_pileup_f func;
57         void *user_data;
58         bam_plbuf_t *plbuf;
59 };
60
61 void bam_lplbuf_reset(bam_lplbuf_t *buf)
62 {
63         freenode_t *p, *q;
64         bam_plbuf_reset(buf->plbuf);
65         for (p = buf->head; p->next;) {
66                 q = p->next;
67                 mp_free(buf->mp, p);
68                 p = q;
69         }
70         buf->head = buf->tail;
71         buf->max_level = 0;
72         buf->n_cur = buf->n_pre = 0;
73         buf->n_nodes = 0;
74 }
75
76 static int tview_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data)
77 {
78         bam_lplbuf_t *tv = (bam_lplbuf_t*)data;
79         freenode_t *p;
80         int i, l, max_level;
81         // allocate memory if necessary
82         if (tv->max < n) { // enlarge
83                 tv->max = n;
84                 kroundup32(tv->max);
85                 tv->cur_level = (int*)realloc(tv->cur_level, sizeof(int) * tv->max);
86                 tv->pre_level = (int*)realloc(tv->pre_level, sizeof(int) * tv->max);
87         }
88         tv->n_cur = n;
89         // update cnt
90         for (p = tv->head; p->next; p = p->next)
91                 if (p->cnt > 0) --p->cnt;
92         // calculate cur_level[]
93         max_level = 0;
94         for (i = l = 0; i < n; ++i) {
95                 const bam_pileup1_t *p = pl + i;
96                 if (p->is_head) {
97                         if (tv->head->next && tv->head->cnt == 0) { // then take a free slot
98                                 freenode_t *p = tv->head->next;
99                                 tv->cur_level[i] = tv->head->level;
100                                 mp_free(tv->mp, tv->head);
101                                 tv->head = p;
102                                 --tv->n_nodes;
103                         } else tv->cur_level[i] = ++tv->max_level;
104                 } else {
105                         tv->cur_level[i] = tv->pre_level[l++];
106                         if (p->is_tail) { // then return a free slot
107                                 tv->tail->level = tv->cur_level[i];
108                                 tv->tail->next = mp_alloc(tv->mp);
109                                 tv->tail = tv->tail->next;
110                                 ++tv->n_nodes;
111                         }
112                 }
113                 if (tv->cur_level[i] > max_level) max_level = tv->cur_level[i];
114                 ((bam_pileup1_t*)p)->level = tv->cur_level[i];
115         }
116         assert(l == tv->n_pre);
117         tv->func(tid, pos, n, pl, tv->user_data);
118         // sort the linked list
119         if (tv->n_nodes) {
120                 freenode_t *q;
121                 if (tv->n_nodes + 1 > tv->m_aux) { // enlarge
122                         tv->m_aux = tv->n_nodes + 1;
123                         kroundup32(tv->m_aux);
124                         tv->aux = (freenode_t**)realloc(tv->aux, sizeof(void*) * tv->m_aux);
125                 }
126                 for (p = tv->head, i = l = 0; p->next;) {
127                         if (p->level > max_level) { // then discard this entry
128                                 q = p->next;
129                                 mp_free(tv->mp, p);
130                                 p = q;
131                         } else {
132                                 tv->aux[i++] = p;
133                                 p = p->next;
134                         }
135                 }
136                 tv->aux[i] = tv->tail; // add a proper tail for the loop below
137                 tv->n_nodes = i;
138                 if (tv->n_nodes) {
139                         ks_introsort(node, tv->n_nodes, tv->aux);
140                         for (i = 0; i < tv->n_nodes; ++i) tv->aux[i]->next = tv->aux[i+1];
141                         tv->head = tv->aux[0];
142                 } else tv->head = tv->tail;
143         }
144         // clean up
145         tv->max_level = max_level;
146         memcpy(tv->pre_level, tv->cur_level, tv->n_cur * 4);
147         // squeeze out terminated levels
148         for (i = l = 0; i < n; ++i) {
149                 const bam_pileup1_t *p = pl + i;
150                 if (!p->is_tail)
151                         tv->pre_level[l++] = tv->pre_level[i];
152         }
153         tv->n_pre = l;
154 /*
155         fprintf(stderr, "%d\t", pos+1);
156         for (i = 0; i < n; ++i) {
157                 const bam_pileup1_t *p = pl + i;
158                 if (p->is_head) fprintf(stderr, "^");
159                 if (p->is_tail) fprintf(stderr, "$");
160                 fprintf(stderr, "%d,", p->level);
161         }
162         fprintf(stderr, "\n");
163 */
164         return 0;
165 }
166
167 bam_lplbuf_t *bam_lplbuf_init(bam_pileup_f func, void *data)
168 {
169         bam_lplbuf_t *tv;
170         tv = (bam_lplbuf_t*)calloc(1, sizeof(bam_lplbuf_t));
171         tv->mp = mp_init();
172         tv->head = tv->tail = mp_alloc(tv->mp);
173         tv->func = func;
174         tv->user_data = data;
175         tv->plbuf = bam_plbuf_init(tview_func, tv);
176         return (bam_lplbuf_t*)tv;
177 }
178
179 void bam_lplbuf_destroy(bam_lplbuf_t *tv)
180 {
181         freenode_t *p, *q;
182         free(tv->cur_level); free(tv->pre_level);
183         bam_plbuf_destroy(tv->plbuf);
184         free(tv->aux);
185         for (p = tv->head; p->next;) {
186                 q = p->next;
187                 mp_free(tv->mp, p); p = q;
188         }
189         mp_free(tv->mp, p);
190         assert(tv->mp->cnt == 0);
191         mp_destroy(tv->mp);
192         free(tv);
193 }
194
195 int bam_lplbuf_push(const bam1_t *b, bam_lplbuf_t *tv)
196 {
197         return bam_plbuf_push(b, tv->plbuf);
198 }