StochHMM  v0.34
Flexible Hidden Markov Model C++ Library and Application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
stoch_viterbi.cpp
Go to the documentation of this file.
1 //
2 // stoch_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  void trellis::stochastic_viterbi(model* h, sequences *sqs){
14  hmm = h;
15  seqs = sqs;
16  seq_size = seqs->getLength();
17  state_size = hmm->state_size();
18  exDef_defined = seqs->exDefDefined();
19 
20 
21  stochastic_viterbi();
22 
23  }
24 
25  void trellis::stochastic_viterbi(){
26  if (hmm->isBasic()){
27  simple_stochastic_viterbi();
28  }
29  }
30 
31 
32  void trellis::simple_stochastic_viterbi(model* h, sequences* sqs){
33  hmm = h;
34  seqs = sqs;
35  seq_size = seqs->getLength();
36  state_size = hmm->state_size();
37  exDef_defined = seqs->exDefDefined();
38 
39  simple_stochastic_viterbi();
40  }
41 
42  void trellis::simple_stochastic_viterbi(){
43  scoring_previous = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
44  scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
45  traceback_table = new int_2D(seq_size,std::vector<int16_t> (state_size,-1));
46  stochastic_table = new (std::nothrow) stochTable(seq_size);
47 
48  std::bitset<STATE_MAX> next_states;
49  std::bitset<STATE_MAX> current_states;
50 
51  double viterbi_temp(-INFINITY);
52  double emission(-INFINITY);
53  bool exDef_position(false);
54  ending_viterbi_score = -INFINITY;
55 
56  state* init = hmm->getInitial();
57 
58  std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
59  std::bitset<STATE_MAX>* from_trans(NULL);
60 
61  //Calculate Viterbi from transitions from INIT (initial) state
62  for(size_t st = 0; st < state_size; ++st){
63  if ((*initial_to)[st]){ //if the bitset is set (meaning there is a transition to this state), calculate the viterbi
64 
65  viterbi_temp = (*hmm)[st]->get_emission_prob(*seqs,0) + getTransition(init, st, 0);;
66 
67  if (viterbi_temp > -INFINITY){
68  if ((*scoring_current)[st] < viterbi_temp){
69  (*scoring_current)[st] = viterbi_temp;
70  }
71  next_states |= (*(*hmm)[st]->getTo());
72  }
73  }
74  }
75 
76  for(size_t position = 1; position < seq_size ; ++position ){
77 
78  //Swap current and previous viterbi scores
79  scoring_previous->assign(state_size,-INFINITY);
80  swap_ptr = scoring_previous;
81  scoring_previous = scoring_current;
82  scoring_current = swap_ptr;
83 
84  //Swap current_states and next states sets
85  current_states.reset();
86  current_states |= next_states;
87  next_states.reset();
88 
89  if (exDef_defined){
90  exDef_position = seqs->exDefDefined(position);
91  }
92 
93  for (size_t st_current = 0; st_current < state_size; ++st_current){ //i is current state that emits value
94 
95  //Check to see if current state is valid
96  if (!current_states[st_current]){
97  continue;
98  }
99 
100  emission = (*hmm)[st_current]->get_emission_prob(*seqs, position);
101 
102 
103  if (exDef_defined && exDef_position){
104  emission += seqs->getWeight(position, st_current);
105  }
106 
107  if (emission == -INFINITY){
108  continue;
109  }
110 
111  //Get list of states that are valid previous states
112  from_trans = (*hmm)[st_current]->getFrom();
113 
114  for (size_t st_previous = 0; st_previous < state_size ; ++st_previous) { //j is previous state
115  if (!(*from_trans)[st_previous]){
116  continue;
117  }
118 
119  if ((*scoring_previous)[st_previous] != -INFINITY){
120  viterbi_temp = getTransition((*hmm)[st_previous], st_current , position) + emission + (*scoring_previous)[st_previous];
121 
122  if (viterbi_temp == -INFINITY){
123  continue;
124  }
125 
126  //Save partial value to stochastic table
127  stochastic_table->push(position-1,st_current,st_previous,viterbi_temp);
128 
129  if (viterbi_temp > (*scoring_current)[st_current]){
130  (*scoring_current)[st_current] = viterbi_temp;
131  (*traceback_table)[position][st_current] = st_previous;
132  }
133 
134  next_states |= (*(*hmm)[st_current]->getTo());
135  }
136  }
137  }
138  }
139 
140  //TODO: Calculate ending and set the final viterbi and traceback pointer
141  //Swap current and previous viterbi scores
142  scoring_previous->assign(state_size,-INFINITY);
143  swap_ptr = scoring_previous;
144  scoring_previous = scoring_current;
145  scoring_current = swap_ptr;
146 
147  for(size_t st_previous = 0; st_previous < state_size ;++st_previous){
148  if ((*scoring_previous)[st_previous] > -INFINITY){
149  viterbi_temp = (*scoring_previous)[st_previous] + (*hmm)[st_previous]->getEndTrans();
150  if (viterbi_temp == -INFINITY){
151  continue;
152  }
153  stochastic_table->push(seq_size-1,SIZE_MAX, st_previous,viterbi_temp);
154 
155  if (viterbi_temp > ending_viterbi_score){
156  ending_viterbi_score = viterbi_temp;
157  ending_viterbi_tb = st_previous;
158  }
159  }
160  }
161 
162  stochastic_table->finalize();
163  //stochastic_table->print();
164 
165  delete scoring_previous;
166  delete scoring_current;
167  scoring_current = NULL;
168  scoring_previous = NULL;
169  }
170 
171  void trellis::naive_stochastic_viterbi(model* h, sequences* sqs){
172  hmm = h;
173  seqs = sqs;
174  seq_size = seqs->getLength();
175  state_size = hmm->state_size();
176  exDef_defined = seqs->exDefDefined();
177 
178  naive_stochastic_viterbi();
179  }
180 
181  void trellis::naive_stochastic_viterbi(){
182  traceback_table = new(std::nothrow) int_2D(seq_size, std::vector<int16_t>(state_size,-1));
183  dbl_viterbi_score = new (std::nothrow) double_2D(seq_size, std::vector<double>(state_size,-INFINITY));
184  stochastic_table = new (std::nothrow) stochTable(seq_size);
185 
186  double emission(-INFINITY);
187  double viterbi_temp(-INFINITY);
188  double trans(-INFINITY);
189  double previous(-INFINITY);
190  bool exDef_position(false);
191  ending_viterbi_tb = -1;
192  ending_viterbi_score = -INFINITY;
193 
194  state* init = hmm->getInitial();
195 
196  //Calculate from Initial states
197  for(size_t st = 0; st < state_size; ++st){
198  viterbi_temp = (*hmm)[st]->get_emission_prob(*seqs,0) + getTransition(init, st, 0);
199  (*dbl_viterbi_score)[0][st]=viterbi_temp;
200  (*traceback_table)[0][st]=-1;
201  }
202 
203  //Calculate Forward for all states
204  for (size_t position = 1 ; position < seq_size ; ++position){
205 
206  if (exDef_defined){
207  exDef_position = seqs->exDefDefined(position);
208  }
209 
210  for (size_t st_current = 0; st_current < state_size; ++st_current){
211  //Calc emissions
212  emission = (*hmm)[st_current]->get_emission_prob(*seqs, position);
213 
214  if (exDef_defined && exDef_position){
215  emission += seqs->getWeight(position, st_current);
216  }
217 
218  if (emission == -INFINITY){
219  continue;
220  }
221 
222  for (size_t st_previous = 0; st_previous < state_size; ++st_previous){
223  previous = (*dbl_viterbi_score)[position-1][st_previous];
224  if (previous == -INFINITY){
225  continue;
226  }
227 
228  trans = getTransition(hmm->getState(st_previous), st_current, position);
229 
230  if (trans !=-INFINITY){
231  viterbi_temp = emission + trans + previous;
232 
233  //Save partial value to stochastic table
234  stochastic_table->push(position-1,st_current,st_previous,viterbi_temp);
235 
236  if (viterbi_temp > (*dbl_viterbi_score)[position][st_current]){
237  (*dbl_viterbi_score)[position][st_current] = viterbi_temp;
238  (*traceback_table)[position][st_current] = previous;
239  }
240  }
241  }
242  }
243  }
244 
245 
246  //Calculate Ending Transition
247  ending_viterbi_tb = -1;
248  for (size_t st_previous = 0; st_previous < state_size; ++st_previous){
249 
250  if ((*hmm)[st_previous]->getEndTrans() != -INFINITY){
251 
252  if ((*dbl_viterbi_score)[seq_size-1][st_previous] != -INFINITY){
253  viterbi_temp = (*dbl_viterbi_score)[seq_size-1][st_previous] + (*hmm)[st_previous]->getEndTrans();
254  stochastic_table->push(seq_size-1,SIZE_MAX, st_previous,viterbi_temp);
255 
256  if (viterbi_temp > ending_viterbi_score){
257  ending_viterbi_score = viterbi_temp;
258  ending_viterbi_tb = st_previous;
259  }
260  }
261  }
262 
263  }
264 
265  stochastic_table->finalize();
266  //stochastic_table->print();
267 
268  return;
269  }
270 
271 // void trellis::simple_alt_stochastic_viterbi(model* h, sequences* sqs){
272 // hmm = h;
273 // seqs = sqs;
274 // seq_size = seqs->getLength();
275 // state_size = hmm->state_size();
276 // exDef_defined = seqs->exDefDefined();
277 //
278 // scoring_previous = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
279 // scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
280 // traceback_table = new int_2D(seq_size,std::vector<int16_t> (state_size,-1));
281 // alt_stochastic_table = new (std::nothrow) alt_stochTable(state_size,seq_size);
282 //
283 // std::bitset<STATE_MAX> next_states;
284 // std::bitset<STATE_MAX> current_states;
285 //
286 // double viterbi_temp(-INFINITY);
287 // double emission(-INFINITY);
288 // bool exDef_position(false);
289 // ending_viterbi_score = -INFINITY;
290 //
291 // state* init = hmm->getInitial();
292 //
293 // std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
294 // std::bitset<STATE_MAX>* from_trans(NULL);
295 //
296 // //Calculate Viterbi from transitions from INIT (initial) state
297 // for(size_t st = 0; st < state_size; ++st){
298 // if ((*initial_to)[st]){ //if the bitset is set (meaning there is a transition to this state), calculate the viterbi
299 //
300 // viterbi_temp = (*hmm)[st]->get_emission_prob(*seqs,0) + getTransition(init, st, 0);;
301 //
302 // if (viterbi_temp > -INFINITY){
303 // if ((*scoring_current)[st] < viterbi_temp){
304 // (*scoring_current)[st] = viterbi_temp;
305 // }
306 // next_states |= (*(*hmm)[st]->getTo());
307 // }
308 // }
309 // }
310 //
311 // for(size_t position = 1; position < seq_size ; ++position ){
312 //
313 // //Swap current and previous viterbi scores
314 // scoring_previous->assign(state_size,-INFINITY);
315 // swap_ptr = scoring_previous;
316 // scoring_previous = scoring_current;
317 // scoring_current = swap_ptr;
318 //
319 // //Swap current_states and next states sets
320 // current_states.reset();
321 // current_states |= next_states;
322 // next_states.reset();
323 //
324 // if (exDef_defined){
325 // exDef_position = seqs->exDefDefined(position);
326 // }
327 //
328 // for (size_t st_current = 0; st_current < state_size; ++st_current){ //i is current state that emits value
329 //
330 // //Check to see if current state is valid
331 // if (!current_states[st_current]){
332 // continue;
333 // }
334 //
335 // emission = (*hmm)[st_current]->get_emission_prob(*seqs, position);
336 //
337 //
338 // if (exDef_defined && exDef_position){
339 // emission += seqs->getWeight(position, st_current);
340 // }
341 //
342 // if (emission == -INFINITY){
343 // continue;
344 // }
345 //
346 // //Get list of states that are valid previous states
347 // from_trans = (*hmm)[st_current]->getFrom();
348 //
349 // for (size_t st_previous = 0; st_previous < state_size ; ++st_previous) { //j is previous state
350 // if (!(*from_trans)[st_previous]){
351 // continue;
352 // }
353 //
354 // if ((*scoring_previous)[st_previous] != -INFINITY){
355 // viterbi_temp = getTransition((*hmm)[st_previous], st_current , position) + emission + (*scoring_previous)[st_previous];
356 //
357 // if (viterbi_temp == -INFINITY){
358 // continue;
359 // }
360 //
361 // //Save partial value to stochastic table
362 // alt_stochastic_table->push(position-1,st_current,st_previous,viterbi_temp);
363 //
364 // if (viterbi_temp > (*scoring_current)[st_current]){
365 // (*scoring_current)[st_current] = viterbi_temp;
366 // (*traceback_table)[position][st_current] = st_previous;
367 // }
368 //
369 // next_states |= (*(*hmm)[st_current]->getTo());
370 // }
371 // }
372 // }
373 // }
374 //
375 // //TODO: Calculate ending and set the final viterbi and traceback pointer
376 // //Swap current and previous viterbi scores
377 // scoring_previous->assign(state_size,-INFINITY);
378 // swap_ptr = scoring_previous;
379 // scoring_previous = scoring_current;
380 // scoring_current = swap_ptr;
381 //
382 // for(size_t st_previous = 0; st_previous < state_size ;++st_previous){
383 // if ((*scoring_previous)[st_previous] > -INFINITY){
384 // viterbi_temp = (*scoring_previous)[st_previous] + (*hmm)[st_previous]->getEndTrans();
385 // if (viterbi_temp == -INFINITY){
386 // continue;
387 // }
388 // alt_stochastic_table->push_ending(st_previous,viterbi_temp);
389 //
390 // if (viterbi_temp > ending_viterbi_score){
391 // ending_viterbi_score = viterbi_temp;
392 // ending_viterbi_tb = st_previous;
393 // }
394 // }
395 // }
396 //
397 // alt_stochastic_table->finalize();
398 // alt_stochastic_table->print();
399 //
400 // delete scoring_previous;
401 // delete scoring_current;
402 // scoring_current = NULL;
403 // scoring_previous = NULL;
404 // }
405 
406 
407  void trellis::simple_simple_stochastic_viterbi(model* h, sequences* sqs){
408  hmm = h;
409  seqs = sqs;
410  seq_size = seqs->getLength();
411  state_size = hmm->state_size();
412  exDef_defined = seqs->exDefDefined();
413 
414  scoring_previous = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
415  scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
416  traceback_table = new int_2D(seq_size,std::vector<int16_t> (state_size,-1));
417  alt_simple_stochastic_table = new (std::nothrow) alt_simple_stochTable(state_size,seq_size);
418 
419  std::bitset<STATE_MAX> next_states;
420  std::bitset<STATE_MAX> current_states;
421 
422  double viterbi_temp(-INFINITY);
423  double emission(-INFINITY);
424  bool exDef_position(false);
425  ending_viterbi_score = -INFINITY;
426 
427  state* init = hmm->getInitial();
428 
429  std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
430  std::bitset<STATE_MAX>* from_trans(NULL);
431 
432  //Calculate Viterbi from transitions from INIT (initial) state
433  for(size_t st = 0; st < state_size; ++st){
434  if ((*initial_to)[st]){ //if the bitset is set (meaning there is a transition to this state), calculate the viterbi
435 
436  viterbi_temp = (*hmm)[st]->get_emission_prob(*seqs,0) + getTransition(init, st, 0);;
437 
438  if (viterbi_temp > -INFINITY){
439  if ((*scoring_current)[st] < viterbi_temp){
440  (*scoring_current)[st] = viterbi_temp;
441  }
442  next_states |= (*(*hmm)[st]->getTo());
443  }
444  }
445  }
446 
447  for(size_t position = 1; position < seq_size ; ++position ){
448 
449  //Swap current and previous viterbi scores
450  scoring_previous->assign(state_size,-INFINITY);
451  swap_ptr = scoring_previous;
452  scoring_previous = scoring_current;
453  scoring_current = swap_ptr;
454 
455  //Swap current_states and next states sets
456  current_states.reset();
457  current_states |= next_states;
458  next_states.reset();
459 
460  if (exDef_defined){
461  exDef_position = seqs->exDefDefined(position);
462  }
463 
464  for (size_t st_current = 0; st_current < state_size; ++st_current){ //i is current state that emits value
465 
466  //Check to see if current state is valid
467  if (!current_states[st_current]){
468  continue;
469  }
470 
471  emission = (*hmm)[st_current]->get_emission_prob(*seqs, position);
472 
473 
474  if (exDef_defined && exDef_position){
475  emission += seqs->getWeight(position, st_current);
476  }
477 
478  if (emission == -INFINITY){
479  continue;
480  }
481 
482  //Get list of states that are valid previous states
483  from_trans = (*hmm)[st_current]->getFrom();
484 
485  for (size_t st_previous = 0; st_previous < state_size ; ++st_previous) { //j is previous state
486  if (!(*from_trans)[st_previous]){
487  continue;
488  }
489 
490  if ((*scoring_previous)[st_previous] != -INFINITY){
491  viterbi_temp = getTransition((*hmm)[st_previous], st_current , position) + emission + (*scoring_previous)[st_previous];
492 
493  if (viterbi_temp == -INFINITY){
494  continue;
495  }
496 
497  //Save partial value to stochastic table
498  alt_simple_stochastic_table->push(position-1,st_current,st_previous,viterbi_temp);
499 
500  if (viterbi_temp > (*scoring_current)[st_current]){
501  (*scoring_current)[st_current] = viterbi_temp;
502  (*traceback_table)[position][st_current] = st_previous;
503  }
504 
505  next_states |= (*(*hmm)[st_current]->getTo());
506  }
507  }
508  }
509  }
510 
511  //TODO: Calculate ending and set the final viterbi and traceback pointer
512  //Swap current and previous viterbi scores
513  scoring_previous->assign(state_size,-INFINITY);
514  swap_ptr = scoring_previous;
515  scoring_previous = scoring_current;
516  scoring_current = swap_ptr;
517 
518  for(size_t st_previous = 0; st_previous < state_size ;++st_previous){
519  if ((*scoring_previous)[st_previous] > -INFINITY){
520  viterbi_temp = (*scoring_previous)[st_previous] + (*hmm)[st_previous]->getEndTrans();
521  if (viterbi_temp == -INFINITY){
522  continue;
523  }
524  alt_simple_stochastic_table->push_ending(st_previous,viterbi_temp);
525 
526  if (viterbi_temp > ending_viterbi_score){
527  ending_viterbi_score = viterbi_temp;
528  ending_viterbi_tb = st_previous;
529  }
530  }
531  }
532 
533  alt_simple_stochastic_table->finalize();
534  //alt_simple_stochastic_table->print();
535 
536  delete scoring_previous;
537  delete scoring_current;
538  scoring_current = NULL;
539  scoring_previous = NULL;
540  }
541 
542 }
Generated on Tue Jul 30 2013 13:23:12 for StochHMM by   doxygen 1.8.1