1 /* The MIT License 2 3 Copyright (c) 2014-2016 Genome Research Ltd. 4 5 Author: Petr Danecek <pd3@sanger.ac.uk> 6 7 Permission is hereby granted, free of charge, to any person obtaining a copy 8 of this software and associated documentation files (the "Software"), to deal 9 in the Software without restriction, including without limitation the rights 10 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 11 copies of the Software, and to permit persons to whom the Software is 12 furnished to do so, subject to the following conditions: 13 14 The above copyright notice and this permission notice shall be included in 15 all copies or substantial portions of the Software. 16 17 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 18 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 19 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 20 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 21 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 22 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 23 THE SOFTWARE. 24 25 */ 26 27 #ifndef __HMM_H__ 28 #define __HMM_H__ 29 30 #define MAT(matrix,ndim,i,j) (matrix)[(ndim)*(i)+(j)] // P(i|j), that is, transition j->i 31 32 typedef struct _hmm_t hmm_t; 33 34 typedef void (*set_tprob_f) (hmm_t *hmm, uint32_t prev_pos, uint32_t pos, void *data, double *tprob); 35 36 /** 37 * hmm_init() - initialize HMM 38 * @nstates: number of states 39 * @tprob: transition probabilities matrix (nstates x nstates), for elements ordering 40 * see the MAT macro above. 41 * @ntprob: number of precalculated tprob matrices or 0 for constant probs, independent 42 * of distance 43 */ 44 hmm_t *hmm_init(int nstates, double *tprob, int ntprob); 45 void hmm_set_tprob(hmm_t *hmm, double *tprob, int ntprob); 46 47 #define HMM_VIT 1 48 #define HMM_FWD 2 49 #define HMM_BWD 4 50 51 /** 52 * hmm_init_states() - initial state probabilities 53 * @probs: initial state probabilities or NULL to reset to default 54 * 55 * If uncalled, all states are initialized with the same likelihood 56 */ 57 void hmm_init_states(hmm_t *hmm, double *probs); 58 59 /** 60 * hmm_snapshot() - take the model's snapshot, intended for sliding HMM 61 * @snapshot: NULL or snapshot returned by previous hmm_snapshot() call, must be free()-ed by the caller 62 * @pos: take the snapshot at this position 63 * 64 * If both restore() and snapshot() are needed, restore() must be called first. 65 */ 66 void *hmm_snapshot(hmm_t *hmm, void *snapshot, uint32_t pos); 67 68 /** 69 * hmm_restore() - restore model's snapshot, intended for sliding HMM 70 * @snapshot: snapshot returned by hmm_snapshot() call or NULL to reset 71 * @isite: take the snapshot at i-th step 72 * 73 * If both restore() and snapshot() are needed, restore() must be called first. 74 */ 75 void hmm_restore(hmm_t *hmm, void *snapshot); 76 void hmm_reset(hmm_t *hmm, void *snapshot); 77 78 /** 79 * hmm_get_tprob() - return the array of transition matrices, precalculated 80 * to ntprob positions. The first matrix is the initial tprob matrix 81 * set by hmm_init() or hmm_set_tprob() 82 */ 83 double *hmm_get_tprob(hmm_t *hmm); 84 int hmm_get_nstates(hmm_t *hmm); 85 86 /** 87 * hmm_set_tprob_func() - custom setter of transition probabilities 88 */ 89 void hmm_set_tprob_func(hmm_t *hmm, set_tprob_f set_tprob, void *data); 90 91 /** 92 * hmm_run_viterbi() - run Viterbi algorithm 93 * @nsites: number of sites 94 * @eprob: emission probabilities for each site and state (nsites x nstates) 95 * @sites: list of positions 96 * 97 * When done, hmm->vpath[] contains the calculated Viterbi path. The states 98 * are indexed starting from 0, a state at i-th site can be accessed as 99 * vpath[nstates*i]. 100 */ 101 void hmm_run_viterbi(hmm_t *hmm, int nsites, double *eprob, uint32_t *sites); 102 103 /** 104 * hmm_get_viterbi_path() - the viterbi path: state at ith site is the 105 * (nstates*isite)-th element 106 */ 107 uint8_t *hmm_get_viterbi_path(hmm_t *hmm); 108 109 /** 110 * hmm_run_fwd_bwd() - run the forward-backward algorithm 111 * @nsites: number of sites 112 * @eprob: emission probabilities for each site and state (nsites x nstates) 113 * @sites: list of positions 114 */ 115 void hmm_run_fwd_bwd(hmm_t *hmm, int nsites, double *eprob, uint32_t *sites); 116 117 /** 118 * hmm_get_fwd_bwd_prob() - the probability of i-th state at j-th site can 119 * be accessed as fwd_bwd[j*nstates+i]. 120 */ 121 double *hmm_get_fwd_bwd_prob(hmm_t *hmm); 122 123 /** 124 * hmm_run_baum_welch() - run one iteration of Baum-Welch algorithm 125 * @nsites: number of sites 126 * @eprob: emission probabilities for each site and state (nsites x nstates) 127 * @sites: list of positions 128 * 129 * Same as hmm_run_fwd_bwd, in addition a pointer to a matrix with the new 130 * transition probabilities is returned. In this verison, emission 131 * probabilities are not updated. 132 */ 133 double *hmm_run_baum_welch(hmm_t *hmm, int nsites, double *eprob, uint32_t *sites); 134 135 void hmm_destroy(hmm_t *hmm); 136 137 #endif 138 139