1 #include "bcftools.pysam.h"
2
3 /* The MIT License
4
5 Copyright (c) 2014-2017 Genome Research Ltd.
6
7 Author: Petr Danecek <pd3@sanger.ac.uk>
8
9 Permission is hereby granted, free of charge, to any person obtaining a copy
10 of this software and associated documentation files (the "Software"), to deal
11 in the Software without restriction, including without limitation the rights
12 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13 copies of the Software, and to permit persons to whom the Software is
14 furnished to do so, subject to the following conditions:
15
16 The above copyright notice and this permission notice shall be included in
17 all copies or substantial portions of the Software.
18
19 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
25 THE SOFTWARE.
26
27 */
28
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <string.h>
32 #include <assert.h>
33 #include <htslib/hts.h>
34 #include "HMM.h"
35
36 typedef struct
37 {
38 int nstates; // number of hmm's states
39 uint32_t snap_at_pos; // snapshot at this position, 0 when inactive
40 double *vit_prob; // viterbi probabilities, NULL for uniform probs
41 double *fwd_prob; // transition probabilities
42 double *bwd_prob; // transition probabilities
43 }
44 snapshot_t;
45
46 struct _hmm_t
47 {
48 int nstates; // number of states
49
50 double *vprob, *vprob_tmp; // viterbi probs [nstates]
51 uint8_t *vpath; // viterbi path [nstates*nvpath]
52 double *bwd, *bwd_tmp; // bwd probs [nstates]
53 double *fwd; // fwd probs [nstates*(nfwd+1)]
54 int nvpath, nfwd;
55
56 int ntprob_arr; // number of pre-calculated tprob matrices
57 double *curr_tprob, *tmp; // Temporary arrays; curr_tprob is short lived, valid only for
58 // one site (that is, one step of Viterbi algorithm)
59 double *tprob_arr; // Array of transition matrices, precalculated to ntprob_arr
60 // positions. The first matrix is the initial tprob matrix
61 // set by hmm_init() or hmm_set_tprob()
62 set_tprob_f set_tprob; // Optional user function to set / modify transition probabilities
63 // at each site (one step of Viterbi algorithm)
64 void *set_tprob_data;
65 snapshot_t init, state; // Initial and current state probs. Set state from snapshot if prev_snap_pos!=0 or from init otherwise
66 snapshot_t *snapshot; // snapshot->snap_at_pos .. request a snapshot at this position
67 // hmm->state.snap_at_pos .. the current state comes from snapshot made at this position
68 };
69
hmm_get_viterbi_path(hmm_t * hmm)70 uint8_t *hmm_get_viterbi_path(hmm_t *hmm) { return hmm->vpath; }
hmm_get_tprob(hmm_t * hmm)71 double *hmm_get_tprob(hmm_t *hmm) { return hmm->tprob_arr; }
hmm_get_nstates(hmm_t * hmm)72 int hmm_get_nstates(hmm_t *hmm) { return hmm->nstates; }
hmm_get_fwd_bwd_prob(hmm_t * hmm)73 double *hmm_get_fwd_bwd_prob(hmm_t *hmm) { return hmm->fwd; }
74
multiply_matrix(int n,double * a,double * b,double * dst,double * tmp)75 static inline void multiply_matrix(int n, double *a, double *b, double *dst, double *tmp)
76 {
77 double *out = dst;
78 if ( a==dst || b==dst )
79 out = tmp;
80
81 int i,j,k;
82 for (i=0; i<n; i++)
83 {
84 for (j=0; j<n; j++)
85 {
86 double val = 0;
87 for (k=0; k<n; k++) val += MAT(a,n,i,k)*MAT(b,n,k,j);
88 MAT(out,n,i,j) = val;
89 }
90 }
91 if ( out!=dst )
92 memcpy(dst,out,sizeof(double)*n*n);
93 }
94
hmm_init_states(hmm_t * hmm,double * probs)95 void hmm_init_states(hmm_t *hmm, double *probs)
96 {
97 hmm->init.snap_at_pos = hmm->state.snap_at_pos = 0;
98
99 if ( !hmm->init.vit_prob )
100 hmm->init.vit_prob = (double*) malloc(sizeof(double)*hmm->nstates);
101 if ( !hmm->init.fwd_prob )
102 hmm->init.fwd_prob = (double*) malloc(sizeof(double)*hmm->nstates);
103 if ( !hmm->init.bwd_prob )
104 hmm->init.bwd_prob = (double*) malloc(sizeof(double)*hmm->nstates);
105
106 if ( !hmm->state.vit_prob )
107 hmm->state.vit_prob = (double*) malloc(sizeof(double)*hmm->nstates);
108 if ( !hmm->state.fwd_prob )
109 hmm->state.fwd_prob = (double*) malloc(sizeof(double)*hmm->nstates);
110 if ( !hmm->state.bwd_prob )
111 hmm->state.bwd_prob = (double*) malloc(sizeof(double)*hmm->nstates);
112
113 int i;
114 if ( probs )
115 {
116 memcpy(hmm->init.vit_prob,probs,sizeof(double)*hmm->nstates);
117 double sum = 0;
118 for (i=0; i<hmm->nstates; i++) sum += hmm->init.vit_prob[i];
119 for (i=0; i<hmm->nstates; i++) hmm->init.vit_prob[i] /= sum;
120 }
121 else
122 for (i=0; i<hmm->nstates; i++) hmm->init.vit_prob[i] = 1./hmm->nstates;
123
124 memcpy(hmm->init.fwd_prob,hmm->init.vit_prob,sizeof(double)*hmm->nstates); // these remain unchanged
125 memcpy(hmm->init.bwd_prob,hmm->init.vit_prob,sizeof(double)*hmm->nstates);
126 memcpy(hmm->state.vit_prob,hmm->init.vit_prob,sizeof(double)*hmm->nstates); // can be changed by snapshotting
127 memcpy(hmm->state.fwd_prob,hmm->init.fwd_prob,sizeof(double)*hmm->nstates);
128 memcpy(hmm->state.bwd_prob,hmm->init.bwd_prob,sizeof(double)*hmm->nstates);
129 }
hmm_init(int nstates,double * tprob,int ntprob)130 hmm_t *hmm_init(int nstates, double *tprob, int ntprob)
131 {
132 hmm_t *hmm = (hmm_t*) calloc(1,sizeof(hmm_t));
133 hmm->nstates = nstates;
134 hmm->curr_tprob = (double*) malloc(sizeof(double)*nstates*nstates);
135 hmm->tmp = (double*) malloc(sizeof(double)*nstates*nstates);
136 hmm_set_tprob(hmm, tprob, ntprob);
137 hmm_init_states(hmm, NULL);
138 return hmm;
139 }
140
hmm_snapshot(hmm_t * hmm,void * _snapshot,uint32_t pos)141 void *hmm_snapshot(hmm_t *hmm, void *_snapshot, uint32_t pos)
142 {
143 snapshot_t *snapshot = (snapshot_t*) _snapshot;
144 if ( snapshot && snapshot->nstates!=hmm->nstates )
145 {
146 free(snapshot);
147 snapshot = NULL;
148 }
149 if ( !snapshot )
150 {
151 // Allocate the snapshot as a single memory block so that it can be
152 // free()-ed by the user. So make sure the arrays are aligned..
153 size_t str_size = sizeof(snapshot_t);
154 size_t dbl_size = sizeof(double);
155 size_t pad_size = (dbl_size - str_size % dbl_size) % dbl_size;
156 uint8_t *mem = (uint8_t*) malloc(str_size + pad_size + dbl_size*2*hmm->nstates);
157 snapshot = (snapshot_t*) mem;
158 snapshot->nstates = hmm->nstates;
159 snapshot->vit_prob = (double*) (mem + str_size + pad_size);
160 snapshot->fwd_prob = snapshot->vit_prob + hmm->nstates;
161 }
162 snapshot->snap_at_pos = pos;
163 hmm->snapshot = snapshot;
164 return snapshot;
165 }
hmm_restore(hmm_t * hmm,void * _snapshot)166 void hmm_restore(hmm_t *hmm, void *_snapshot)
167 {
168 snapshot_t *snapshot = (snapshot_t*) _snapshot;
169 if ( !snapshot || !snapshot->snap_at_pos )
170 {
171 hmm->state.snap_at_pos = 0;
172 memcpy(hmm->state.vit_prob,hmm->init.vit_prob,sizeof(double)*hmm->nstates);
173 memcpy(hmm->state.fwd_prob,hmm->init.fwd_prob,sizeof(double)*hmm->nstates);
174 }
175 else
176 {
177 hmm->state.snap_at_pos = snapshot->snap_at_pos;
178 memcpy(hmm->state.vit_prob,snapshot->vit_prob,sizeof(double)*hmm->nstates);
179 memcpy(hmm->state.fwd_prob,snapshot->fwd_prob,sizeof(double)*hmm->nstates);
180 }
181 }
hmm_reset(hmm_t * hmm,void * _snapshot)182 void hmm_reset(hmm_t *hmm, void *_snapshot)
183 {
184 snapshot_t *snapshot = (snapshot_t*) _snapshot;
185 if ( snapshot ) snapshot->snap_at_pos = 0;
186 hmm->state.snap_at_pos = 0;
187 memcpy(hmm->state.vit_prob,hmm->init.vit_prob,sizeof(double)*hmm->nstates);
188 memcpy(hmm->state.fwd_prob,hmm->init.fwd_prob,sizeof(double)*hmm->nstates);
189 }
190
hmm_set_tprob(hmm_t * hmm,double * tprob,int ntprob)191 void hmm_set_tprob(hmm_t *hmm, double *tprob, int ntprob)
192 {
193 hmm->ntprob_arr = ntprob;
194 if ( ntprob<=0 ) ntprob = 1;
195
196 if ( !hmm->tprob_arr )
197 hmm->tprob_arr = (double*) malloc(sizeof(double)*hmm->nstates*hmm->nstates*ntprob);
198
199 memcpy(hmm->tprob_arr,tprob,sizeof(double)*hmm->nstates*hmm->nstates);
200
201 int i;
202 for (i=1; i<ntprob; i++)
203 multiply_matrix(hmm->nstates, hmm->tprob_arr, hmm->tprob_arr+(i-1)*hmm->nstates*hmm->nstates, hmm->tprob_arr+i*hmm->nstates*hmm->nstates, hmm->tmp);
204 }
205
hmm_set_tprob_func(hmm_t * hmm,set_tprob_f set_tprob,void * data)206 void hmm_set_tprob_func(hmm_t *hmm, set_tprob_f set_tprob, void *data)
207 {
208 hmm->set_tprob = set_tprob;
209 hmm->set_tprob_data = data;
210 }
211
_set_tprob(hmm_t * hmm,int pos_diff)212 static void _set_tprob(hmm_t *hmm, int pos_diff)
213 {
214 assert( pos_diff>=0 );
215
216 int i, n;
217
218 n = hmm->ntprob_arr ? pos_diff % hmm->ntprob_arr : 0; // n-th precalculated matrix
219 memcpy(hmm->curr_tprob, hmm->tprob_arr+n*hmm->nstates*hmm->nstates, sizeof(*hmm->curr_tprob)*hmm->nstates*hmm->nstates);
220
221 if ( hmm->ntprob_arr > 0 )
222 {
223 n = pos_diff / hmm->ntprob_arr; // number of full blocks to jump
224 for (i=0; i<n; i++)
225 multiply_matrix(hmm->nstates, hmm->tprob_arr+(hmm->ntprob_arr-1)*hmm->nstates*hmm->nstates, hmm->curr_tprob, hmm->curr_tprob, hmm->tmp);
226 }
227 }
228
hmm_run_viterbi(hmm_t * hmm,int n,double * eprobs,uint32_t * sites)229 void hmm_run_viterbi(hmm_t *hmm, int n, double *eprobs, uint32_t *sites)
230 {
231 // Init arrays when run for the first time
232 if ( hmm->nvpath < n )
233 {
234 hmm->nvpath = n;
235 hmm->vpath = (uint8_t*) realloc(hmm->vpath, sizeof(uint8_t)*hmm->nvpath*hmm->nstates);
236 }
237 if ( !hmm->vprob )
238 {
239 hmm->vprob = (double*) malloc(sizeof(double)*hmm->nstates);
240 hmm->vprob_tmp = (double*) malloc(sizeof(double)*hmm->nstates);
241 }
242
243 // Init all states with equal likelihood
244 int i,j, nstates = hmm->nstates;
245 memcpy(hmm->vprob, hmm->state.vit_prob, sizeof(*hmm->state.vit_prob)*nstates);
246 uint32_t prev_pos = hmm->state.snap_at_pos ? hmm->state.snap_at_pos : sites[0];
247
248 // Run Viterbi
249 for (i=0; i<n; i++)
250 {
251 uint8_t *vpath = &hmm->vpath[i*nstates];
252 double *eprob = &eprobs[i*nstates];
253
254 int pos_diff = sites[i] == prev_pos ? 0 : sites[i] - prev_pos - 1;
255 _set_tprob(hmm, pos_diff);
256 if ( hmm->set_tprob ) hmm->set_tprob(hmm, prev_pos, sites[i], hmm->set_tprob_data, hmm->curr_tprob);
257 prev_pos = sites[i];
258
259 double vnorm = 0;
260 for (j=0; j<nstates; j++)
261 {
262 double vmax = 0;
263 int k, k_vmax = 0;
264 for (k=0; k<nstates; k++)
265 {
266 double pval = hmm->vprob[k] * MAT(hmm->curr_tprob,hmm->nstates,j,k);
267 if ( vmax < pval ) { vmax = pval; k_vmax = k; }
268 }
269 vpath[j] = k_vmax;
270 hmm->vprob_tmp[j] = vmax * eprob[j];
271 vnorm += hmm->vprob_tmp[j];
272 }
273 for (j=0; j<nstates; j++) hmm->vprob_tmp[j] /= vnorm;
274 double *tmp = hmm->vprob; hmm->vprob = hmm->vprob_tmp; hmm->vprob_tmp = tmp;
275
276 if ( hmm->snapshot && sites[i]==hmm->snapshot->snap_at_pos )
277 memcpy(hmm->snapshot->vit_prob, hmm->vprob, sizeof(*hmm->vprob)*nstates);
278 }
279
280 // Find the most likely state
281 int iptr = 0;
282 for (i=1; i<nstates; i++)
283 if ( hmm->vprob[iptr] < hmm->vprob[i] ) iptr = i;
284
285 // Trace back the Viterbi path, we are reusing vpath for storing the states (vpath[i*nstates])
286 for (i=n-1; i>=0; i--)
287 {
288 assert( iptr<nstates && hmm->vpath[i*nstates + iptr]<nstates );
289 iptr = hmm->vpath[i*nstates + iptr];
290 hmm->vpath[i*nstates] = iptr; // reusing the array for different purpose here
291 }
292 }
293
hmm_run_fwd_bwd(hmm_t * hmm,int n,double * eprobs,uint32_t * sites)294 void hmm_run_fwd_bwd(hmm_t *hmm, int n, double *eprobs, uint32_t *sites)
295 {
296 // Init arrays when run for the first time
297 if ( hmm->nfwd < n )
298 {
299 hmm->nfwd = n;
300 hmm->fwd = (double*) realloc(hmm->fwd, sizeof(double)*(hmm->nfwd+1)*hmm->nstates);
301 }
302 if ( !hmm->bwd )
303 {
304 hmm->bwd = (double*) malloc(sizeof(double)*hmm->nstates);
305 hmm->bwd_tmp = (double*) malloc(sizeof(double)*hmm->nstates);
306 }
307
308
309 int i,j,k, nstates = hmm->nstates;
310 memcpy(hmm->fwd, hmm->state.fwd_prob, sizeof(*hmm->state.fwd_prob)*nstates);
311 memcpy(hmm->bwd, hmm->state.bwd_prob, sizeof(*hmm->state.bwd_prob)*nstates);
312 uint32_t prev_pos = hmm->state.snap_at_pos ? hmm->state.snap_at_pos : sites[0];
313
314 // Run fwd
315 for (i=0; i<n; i++)
316 {
317 double *fwd_prev = &hmm->fwd[i*nstates];
318 double *fwd = &hmm->fwd[(i+1)*nstates];
319 double *eprob = &eprobs[i*nstates];
320
321 int pos_diff = sites[i] == prev_pos ? 0 : sites[i] - prev_pos - 1;
322
323 _set_tprob(hmm, pos_diff);
324 if ( hmm->set_tprob ) hmm->set_tprob(hmm, prev_pos, sites[i], hmm->set_tprob_data, hmm->curr_tprob);
325 prev_pos = sites[i];
326
327 double norm = 0;
328 for (j=0; j<nstates; j++)
329 {
330 double pval = 0;
331 for (k=0; k<nstates; k++)
332 pval += fwd_prev[k] * MAT(hmm->curr_tprob,hmm->nstates,j,k);
333 fwd[j] = pval * eprob[j];
334 norm += fwd[j];
335 }
336 for (j=0; j<nstates; j++) fwd[j] /= norm;
337
338 if ( hmm->snapshot && sites[i]==hmm->snapshot->snap_at_pos )
339 memcpy(hmm->snapshot->fwd_prob, fwd, sizeof(*fwd)*nstates);
340 }
341
342 // Run bwd
343 double *bwd = hmm->bwd, *bwd_tmp = hmm->bwd_tmp;
344 prev_pos = sites[n-1];
345 for (i=0; i<n; i++)
346 {
347 double *fwd = &hmm->fwd[(n-i)*nstates];
348 double *eprob = &eprobs[(n-i-1)*nstates];
349
350 int pos_diff = sites[n-i-1] == prev_pos ? 0 : prev_pos - sites[n-i-1] - 1;
351
352 _set_tprob(hmm, pos_diff);
353 if ( hmm->set_tprob ) hmm->set_tprob(hmm, sites[n-i-1], prev_pos, hmm->set_tprob_data, hmm->curr_tprob);
354 prev_pos = sites[n-i-1];
355
356 double bwd_norm = 0;
357 for (j=0; j<nstates; j++)
358 {
359 double pval = 0;
360 for (k=0; k<nstates; k++)
361 pval += bwd[k] * eprob[k] * MAT(hmm->curr_tprob,hmm->nstates,k,j);
362 bwd_tmp[j] = pval;
363 bwd_norm += pval;
364 }
365 double norm = 0;
366 for (j=0; j<nstates; j++)
367 {
368 bwd_tmp[j] /= bwd_norm;
369 fwd[j] *= bwd_tmp[j]; // fwd now stores fwd*bwd
370 norm += fwd[j];
371 }
372 for (j=0; j<nstates; j++) fwd[j] /= norm;
373 double *tmp = bwd_tmp; bwd_tmp = bwd; bwd = tmp;
374 }
375 }
376
hmm_run_baum_welch(hmm_t * hmm,int n,double * eprobs,uint32_t * sites)377 double *hmm_run_baum_welch(hmm_t *hmm, int n, double *eprobs, uint32_t *sites)
378 {
379 // Init arrays when run for the first time
380 if ( hmm->nfwd < n )
381 {
382 hmm->nfwd = n;
383 hmm->fwd = (double*) realloc(hmm->fwd, sizeof(double)*(hmm->nfwd+1)*hmm->nstates);
384 }
385 if ( !hmm->bwd )
386 {
387 hmm->bwd = (double*) malloc(sizeof(double)*hmm->nstates);
388 hmm->bwd_tmp = (double*) malloc(sizeof(double)*hmm->nstates);
389 }
390
391 // Init all states with equal likelihood
392 int i,j,k, nstates = hmm->nstates;
393 memcpy(hmm->fwd, hmm->state.fwd_prob, sizeof(*hmm->state.fwd_prob)*nstates);
394 memcpy(hmm->bwd, hmm->state.bwd_prob, sizeof(*hmm->state.bwd_prob)*nstates);
395 uint32_t prev_pos = hmm->state.snap_at_pos ? hmm->state.snap_at_pos : sites[0];
396
397 // New transition matrix: temporary values
398 double *tmp_xi = (double*) calloc(nstates*nstates,sizeof(double));
399 double *tmp_gamma = (double*) calloc(nstates,sizeof(double));
400 double *fwd_bwd = (double*) malloc(sizeof(double)*nstates);
401
402 // Run fwd
403 for (i=0; i<n; i++)
404 {
405 double *fwd_prev = &hmm->fwd[i*nstates];
406 double *fwd = &hmm->fwd[(i+1)*nstates];
407 double *eprob = &eprobs[i*nstates];
408
409 int pos_diff = sites[i] == prev_pos ? 0 : sites[i] - prev_pos - 1;
410
411 _set_tprob(hmm, pos_diff);
412 if ( hmm->set_tprob ) hmm->set_tprob(hmm, prev_pos, sites[i], hmm->set_tprob_data, hmm->curr_tprob);
413 prev_pos = sites[i];
414
415 double norm = 0;
416 for (j=0; j<nstates; j++)
417 {
418 double pval = 0;
419 for (k=0; k<nstates; k++)
420 pval += fwd_prev[k] * MAT(hmm->curr_tprob,hmm->nstates,j,k);
421 fwd[j] = pval * eprob[j];
422 norm += fwd[j];
423 }
424 for (j=0; j<nstates; j++) fwd[j] /= norm;
425 }
426
427 // Run bwd
428 double *bwd = hmm->bwd, *bwd_tmp = hmm->bwd_tmp;
429 prev_pos = sites[n-1];
430 for (i=0; i<n; i++)
431 {
432 double *fwd = &hmm->fwd[(n-i)*nstates];
433 double *eprob = &eprobs[(n-i-1)*nstates];
434
435 int pos_diff = sites[n-i-1] == prev_pos ? 0 : prev_pos - sites[n-i-1] - 1;
436
437 _set_tprob(hmm, pos_diff);
438 if ( hmm->set_tprob ) hmm->set_tprob(hmm, sites[n-i-1], prev_pos, hmm->set_tprob_data, hmm->curr_tprob);
439 prev_pos = sites[n-i-1];
440
441 double bwd_norm = 0;
442 for (j=0; j<nstates; j++)
443 {
444 double pval = 0;
445 for (k=0; k<nstates; k++)
446 pval += bwd[k] * eprob[k] * MAT(hmm->curr_tprob,hmm->nstates,k,j);
447 bwd_tmp[j] = pval;
448 bwd_norm += pval;
449 }
450 double norm = 0;
451 for (j=0; j<nstates; j++)
452 {
453 bwd_tmp[j] /= bwd_norm;
454 fwd_bwd[j] = fwd[j]*bwd_tmp[j];
455 norm += fwd_bwd[j];
456 }
457 for (j=0; j<nstates; j++)
458 {
459 fwd_bwd[j] /= norm;
460 tmp_gamma[j] += fwd_bwd[j];
461 }
462
463 for (j=0; j<nstates; j++)
464 {
465 for (k=0; k<nstates; k++)
466 {
467 MAT(tmp_xi,nstates,k,j) += fwd[j]*bwd[k]*MAT(hmm->tprob_arr,hmm->nstates,k,j)*eprob[k] / norm;
468 }
469 }
470
471 for (j=0; j<nstates; j++) fwd[j] = fwd_bwd[j]; // fwd now stores fwd*bwd
472
473 double *tmp = bwd_tmp; bwd_tmp = bwd; bwd = tmp;
474 }
475 for (j=0; j<nstates; j++)
476 {
477 double norm = 0;
478 for (k=0; k<nstates; k++)
479 {
480 MAT(hmm->curr_tprob,nstates,k,j) = MAT(tmp_xi,nstates,k,j) / tmp_gamma[j];
481 norm += MAT(hmm->curr_tprob,nstates,k,j);
482 }
483 for (k=0; k<nstates; k++)
484 MAT(hmm->curr_tprob,nstates,k,j) /= norm;
485 }
486 free(tmp_gamma);
487 free(tmp_xi);
488 free(fwd_bwd);
489 return hmm->curr_tprob;
490 }
491
hmm_destroy(hmm_t * hmm)492 void hmm_destroy(hmm_t *hmm)
493 {
494 free(hmm->init.vit_prob);
495 free(hmm->init.fwd_prob);
496 free(hmm->init.bwd_prob);
497 free(hmm->state.vit_prob);
498 free(hmm->state.fwd_prob);
499 free(hmm->state.bwd_prob);
500 free(hmm->vprob);
501 free(hmm->vprob_tmp);
502 free(hmm->vpath);
503 free(hmm->curr_tprob);
504 free(hmm->tmp);
505 free(hmm->tprob_arr);
506 free(hmm->fwd);
507 free(hmm->bwd);
508 free(hmm->bwd_tmp);
509 free(hmm);
510 }
511
512