StochHMM
v0.34
Flexible Hidden Markov Model C++ Library and Application
Main Page
Modules
Namespaces
Classes
Files
File List
File Members
All
Classes
Namespaces
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Friends
Macros
Groups
Pages
forward_viterbi.cpp
Go to the documentation of this file.
1
//
2
// forward_viterbi.cpp
3
// StochHMM
4
//
5
// Created by Paul Lott on 2/4/13.
6
// Copyright (c) 2013 Korf Lab, Genome Center, UC Davis, Davis, CA. All rights reserved.
7
//
8
9
#include "
trellis.h
"
10
11
namespace
StochHMM {
12
13
//
14
// void trellis::forward_viterbi(StochHMM::model* h , StochHMM::sequences *sqs){
15
// //Initialize the table
16
// hmm = h;
17
// seqs = sqs;
18
// seq_size = seqs->getLength();
19
// state_size = hmm->state_size();
20
// exDef_defined = seqs->exDefDefined();
21
// }
22
//
23
// void trellis::forward_viterbi(){
24
// if (hmm->isBasic()){
25
// simple_forward_viterbi();
26
// }
27
// else{
28
// fast_complex_forward_viterbi();
29
// }
30
// }
31
//
32
//
33
//
34
//
35
// //!
36
// void trellis::simple_forward_viterbi(){
37
//
38
// if (!hmm->isBasic()){
39
// std::cerr << "Model isn't a simple/basic HMM. Use complex algorithms\n";
40
// return;
41
// }
42
//
43
// //Initialize the traceback table
44
// if (traceback_table != NULL || forward_score != NULL){
45
// delete traceback_table;
46
// delete forward_score;
47
// traceback_table = NULL;
48
// forward_score = NULL;
49
// }
50
//
51
// traceback_table = new int_2D(seq_size,std::vector<int16_t> (state_size,-1));
52
// forward_score = new (std::nothrow) float_2D(seq_size, std::vector<float>(state_size,-INFINITY));
53
// scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
54
// scoring_previous= new (std::nothrow) std::vector<double> (state_size,-INFINITY);
55
// alt_scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
56
// alt_scoring_previous= new (std::nothrow) std::vector<double> (state_size,-INFINITY);
57
//
58
// if (scoring_current == NULL || scoring_previous == NULL ||
59
// alt_scoring_current == NULL || alt_scoring_previous == NULL ||
60
// traceback_table == NULL || forward_score == NULL){
61
// std::cerr << "Can't allocate forward score table. OUT OF MEMORY" << std::endl;
62
// exit(2);
63
// }
64
//
65
// std::bitset<STATE_MAX> next_states;
66
// std::bitset<STATE_MAX> current_states;
67
//
68
// double forward_temp(-INFINITY);
69
// double viterbi_temp(-INFINITY);
70
// double emission(-INFINITY);
71
// bool exDef_position(false);
72
//
73
// state* init = hmm->getInitial();
74
//
75
// std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
76
// std::bitset<STATE_MAX>* from_trans(NULL);
77
//
78
//
79
// //Calculate Viterbi and Forward from transitions from INIT (initial) state
80
// for(size_t i = 0; i < state_size; ++i){
81
// if ((*initial_to)[i]){ //if the bitset is set (meaning there is a transition to this state), calculate the viterbi
82
//
83
// forward_temp = (*hmm)[i]->get_emission_prob(*seqs,0) + getTransition(init, i, 0);
84
//
85
//
86
// if (forward_temp > -INFINITY){
87
//
88
// (*forward_score)[0][i] = forward_temp;
89
// (*alt_scoring_current)[i] = forward_temp;
90
// (*scoring_current)[i] = forward_temp;
91
// next_states |= (*(*hmm)[i]->getTo());
92
// }
93
// }
94
// }
95
//
96
//
97
// for(size_t position = 1; position < seq_size ; ++position ){
98
//
99
//
100
// //Swap current and previous viterbi scores
101
// scoring_previous->assign(state_size,-INFINITY);
102
// swap_ptr = scoring_previous;
103
// scoring_previous = scoring_current;
104
// scoring_current = swap_ptr;
105
//
106
// //Swap current and previous viterbi scores
107
// alt_scoring_previous->assign(state_size,-INFINITY);
108
// swap_ptr = alt_scoring_previous;
109
// alt_scoring_previous = alt_scoring_current;
110
// alt_scoring_current = swap_ptr;
111
//
112
// //Swap current_states and next states sets
113
// current_states.reset();
114
// current_states |= next_states;
115
// next_states.reset();
116
//
117
// if (exDef_defined){
118
// exDef_position = seqs->exDefDefined(position);
119
// }
120
//
121
// //Current State
122
// for (size_t current = 0; current < state_size; ++current){ //current state that emits value
123
//
124
// if (!current_states[current]){
125
// continue;
126
// }
127
//
128
// emission = (*hmm)[current]->get_emission_prob(*seqs, position);
129
//
130
// if (exDef_defined && exDef_position){
131
// emission += seqs->getWeight(position, current);
132
// }
133
//
134
// from_trans = (*hmm)[current]->getFrom();
135
//
136
// for (size_t previous = 0; previous < state_size ; ++previous) { //j is previous state
137
//
138
// //Check that previous state has transition to current state
139
// //and that the previous viterbi score is not -INFINITY
140
// if ((*from_trans)[previous] && (*scoring_previous)[previous] != -INFINITY){
141
//
142
// viterbi_temp = getTransition((*hmm)[previous], current , position) + emission;
143
//
144
// forward_temp = viterbi_temp + (*alt_scoring_previous)[previous];
145
// viterbi_temp += (*scoring_previous)[previous];
146
//
147
// if (viterbi_temp > (*scoring_current)[current]){
148
// (*scoring_current)[current] = viterbi_temp;
149
// (*traceback_table)[position][current] = previous;
150
// }
151
//
152
// if ((*alt_scoring_current)[current] == -INFINITY){
153
// (*alt_scoring_current)[current] = forward_temp;
154
// (*forward_score)[position][current] = forward_temp;
155
// }
156
// else{
157
// (*alt_scoring_current)[current] = addLog(forward_temp, (*alt_scoring_current)[current]);
158
// (*forward_score)[position][current] = (*alt_scoring_current)[current]; }
159
//
160
// next_states |= (*(*hmm)[current]->getTo());
161
// }
162
// }
163
// // std::cout << "State: " << current <<"\t" << exp((*forward_score)[position][current]) << std::endl;
164
// }
165
// }
166
//
167
// //Swap current and previous scores
168
// scoring_previous->assign(state_size,-INFINITY);
169
// swap_ptr = scoring_previous;
170
// scoring_previous = scoring_current;
171
// scoring_current = swap_ptr;
172
//
173
// //Swap current and previous scores
174
// alt_scoring_previous->assign(state_size,-INFINITY);
175
// swap_ptr = alt_scoring_previous;
176
// alt_scoring_previous = alt_scoring_current;
177
// alt_scoring_current = swap_ptr;
178
//
179
//
180
// ending_forward_prob = -INFINITY;
181
// for(size_t i = 0; i < state_size ;++i){
182
// //Ending Forward Transition
183
// if ((*alt_scoring_previous)[i] != -INFINITY){
184
// forward_temp = (*alt_scoring_previous)[i] + (*hmm)[i]->getEndTrans();
185
//
186
// if (forward_temp > -INFINITY){
187
// if (ending_forward_prob == -INFINITY){
188
// ending_forward_prob = forward_temp;
189
// }
190
// else{
191
// ending_forward_prob = addLog(ending_forward_prob,forward_temp);
192
// }
193
// }
194
// }
195
//
196
//
197
// //Ending Viterbi Transition
198
// if ((*scoring_previous)[i] > -INFINITY){
199
// viterbi_temp = (*scoring_previous)[i] + (*hmm)[i]->getEndTrans();
200
// if (viterbi_temp > ending_viterbi_score){
201
// ending_viterbi_score = viterbi_temp;
202
// ending_viterbi_tb = i;
203
// }
204
// }
205
//
206
// }
207
//
208
// // std::cout << exp(ending_posterior) << std::endl;
209
// delete scoring_previous;
210
// delete scoring_current;
211
// delete alt_scoring_current;
212
// delete alt_scoring_previous;
213
// scoring_previous = NULL;
214
// scoring_current = NULL;
215
// alt_scoring_current = NULL;
216
// alt_scoring_previous = NULL;
217
// }
218
219
220
221
// void trellis::forward_viterbi(){
222
//
223
// //Initialize the traceback table
224
// if (traceback_table != NULL){
225
// delete traceback_table;
226
// }
227
//
228
// traceback_table = new int_2D(seq_size,std::vector<int16_t> (state_size,1));
229
//
230
// //Storing Viterbi Scores
231
// scoring_previous = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
232
// scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
233
//
234
// //Storing Forward Scores
235
// alt_scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
236
// alt_scoring_previous = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
237
//
238
//
239
// if (scoring_previous == NULL || scoring_current == NULL || alt_scoring_current == NULL ||
240
// alt_scoring_previous == NULL || traceback_table == NULL){
241
// std::cerr << "Can't allocate Viterbi score and traceback table. OUT OF MEMORY" << std::endl;
242
// exit(2);
243
// }
244
//
245
// std::bitset<STATE_MAX> next_states;
246
// std::bitset<STATE_MAX> current_states;
247
//
248
//
249
//
250
// double forward_temp(-INFINITY);
251
// double emission(-INFINITY);
252
// bool exDef_position(false);
253
//
254
// state* init = hmm->getInitial();
255
// std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
256
// std::bitset<STATE_MAX>* from_trans(NULL);
257
//
258
//
259
// // std::cout << "Position: 0" << std::endl;
260
// //Calculate Viterbi from transitions from INIT (initial) state
261
// for(size_t i = 0; i < state_size; ++i){
262
// if ((*initial_to)[i]){ //if the bitset is set (meaning there is a transition to this state), calculate the viterbi
263
//
264
// forward_temp = (*hmm)[i]->get_emission_prob(*seqs,0) + getTransition(init, i, 0);
265
//
266
// if (forward_temp > -INFINITY){
267
//
268
// (*forward_score)[0][i] = forward_temp;
269
// // std::cout << "State: " << i << "\t" << exp(forward_temp) << std::endl;
270
// next_states |= (*(*hmm)[i]->getTo());
271
// }
272
// }
273
// }
274
//
275
// for(size_t position = 1; position < seq_size ; ++position ){
276
//
277
// //Swap current_states and next states sets
278
//
279
// current_states.reset();
280
// current_states |= next_states;
281
// next_states.reset();
282
//
283
// if (exDef_defined){
284
// exDef_position = seqs->exDefDefined(position);
285
// }
286
//
287
// // std::cout << "\nPosition: " << position << std::endl;
288
//
289
// for (size_t i = 0; i < state_size; ++i){ //i is current state that emits value
290
// if (!current_states[i]){
291
// continue;
292
// }
293
//
294
// emission = (*hmm)[i]->get_emission_prob(*seqs, position);
295
//
296
// if (exDef_defined && exDef_position){
297
// emission += seqs->getWeight(position, i);
298
// }
299
//
300
// from_trans = (*hmm)[i]->getFrom();
301
//
302
// for (size_t j = 0; j < state_size ; ++j) { //j is previous state
303
// if (!(*from_trans)[j]){
304
// continue;
305
// }
306
//
307
// if ((*forward_score)[position-1][j] != INFINITY){
308
// forward_temp = getTransition((*hmm)[j], i , position) + emission + (*forward_score)[position-1][j];
309
//
310
// if ((*forward_score)[position][i] == -INFINITY){
311
// (*forward_score)[position][i] = forward_temp;
312
// }
313
// else{
314
// (*forward_score)[position][i] = addLog((double)forward_temp, (double)(*forward_score)[position][i]);
315
// }
316
//
317
// next_states |= (*(*hmm)[i]->getTo());
318
// }
319
// }
320
// // std::cout << "State: " << i <<"\t" << exp((*forward_score)[position][i]) << std::endl;
321
// }
322
//
323
// }
324
//
325
// for(size_t i = 0; i < state_size ;++i){
326
// if ((*forward_score)[seq_size-1][i] > -INFINITY){
327
// forward_temp = (*forward_score)[seq_size-1][i] + (*hmm)[i]->getEndTrans();
328
//
329
// if (forward_temp > -INFINITY){
330
// if (ending_forward_prob == -INFINITY){
331
// ending_forward_prob = forward_temp;
332
// }
333
// else{
334
// ending_forward_prob = addLog(ending_forward_prob,forward_temp);
335
// }
336
// }
337
// }
338
// }
339
// // std::cout << exp(ending_posterior) << std::endl;
340
// }
341
342
343
}
Generated on Tue Jul 30 2013 13:23:11 for StochHMM by
1.8.1
1.8.1
