StochHMM  v0.34
Flexible Hidden Markov Model C++ Library and Application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
viterbi.cpp
Go to the documentation of this file.
1 //
2 // 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::viterbi(model* h, sequences* sqs){
14  //Initialize the table
15  hmm = h;
16  seqs = sqs;
17  seq_size = seqs->getLength();
18  state_size = hmm->state_size();
19  exDef_defined = seqs->exDefDefined();
20 
21  //TODO: determine which model and chose the type of algorithm to use;
22  viterbi();
23  }
24 
25  void trellis::viterbi(){
26  if (hmm->isBasic()){
27  simple_viterbi();
28  }
29  else{
30  fast_complex_viterbi();
31  }
32  }
33 
34  void trellis::simple_viterbi(model* h, sequences* sqs){
35  //Initialize the table
36  hmm = h;
37  seqs = sqs;
38  seq_size = seqs->getLength();
39  state_size = hmm->state_size();
40  exDef_defined = seqs->exDefDefined();
41 
42  simple_viterbi();
43  }
44 
45  void trellis::simple_viterbi(){
46 
47  if (!hmm->isBasic()){
48  std::cerr << "Model isn't a simple/basic HMM. Use complex algorithms\n";
49  return;
50  }
51 
52  //Initialize the traceback table
53  if (traceback_table != NULL){
54  delete traceback_table;
55  }
56 
57  traceback_table = new int_2D(seq_size,std::vector<int16_t> (state_size,-1));
58  scoring_previous = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
59  scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
60 
61  if (scoring_previous == NULL || scoring_current == NULL || traceback_table == NULL){
62  std::cerr << "Can't allocate Viterbi score and traceback table. OUT OF MEMORY" << std::endl;
63  exit(2);
64  }
65 
66 
67  std::bitset<STATE_MAX> next_states;
68  std::bitset<STATE_MAX> current_states;
69 
70  double viterbi_temp(-INFINITY);
71  double emission(-INFINITY);
72  bool exDef_position(false);
73  ending_viterbi_tb = -1;
74  ending_viterbi_score = -INFINITY;
75 
76  state* init = hmm->getInitial();
77 
78  std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
79  std::bitset<STATE_MAX>* from_trans(NULL);
80 
81  //Calculate Viterbi from transitions from INIT (initial) state
82  for(size_t st = 0; st < state_size; ++st){
83  if ((*initial_to)[st]){ //if the bitset is set (meaning there is a transition to this state), calculate the viterbi
84 
85  viterbi_temp = (*hmm)[st]->get_emission_prob(*seqs,0) + getTransition(init, st, 0);
86 
87  if (viterbi_temp > -INFINITY){
88  if ((*scoring_current)[st] < viterbi_temp){
89  (*scoring_current)[st] = viterbi_temp;
90  }
91  next_states |= (*(*hmm)[st]->getTo());
92  }
93  }
94  }
95 
96  //Each position in the sequence
97  for(size_t position = 1; position < seq_size ; ++position ){
98 
99  //Swap current and previous viterbi scores
100  scoring_previous->assign(state_size,-INFINITY);
101  swap_ptr = scoring_previous;
102  scoring_previous = scoring_current;
103  scoring_current = swap_ptr;
104 
105  //Swap current_states and next states sets
106 
107  current_states.reset();
108  current_states |= next_states;
109  next_states.reset();
110 
111  if (exDef_defined){
112  exDef_position = seqs->exDefDefined(position);
113  }
114 
115  //Current states
116  for (size_t st_current = 0; st_current < state_size; ++st_current){ //Current state that emits value
117 
118  //Check to see if current state is valid
119  if (!current_states[st_current]){
120  continue;
121  }
122 
123  //Get emission of current state
124  emission = (*hmm)[st_current]->get_emission_prob(*seqs, position);
125 
126 
127  if (exDef_defined && exDef_position){
128  emission += seqs->getWeight(position, st_current);
129  }
130 
131 
132  if (emission == -INFINITY){
133  continue;
134  }
135 
136  //Get list of states that are valid previous states
137  from_trans = (*hmm)[st_current]->getFrom();
138 
139  for (size_t st_previous = 0; st_previous < state_size ; ++st_previous) { //for previous states
140  if (!(*from_trans)[st_previous]){
141  continue;
142  }
143 
144  //Check that previous state has transition to current state
145  //and that the previous viterbi score is not -INFINITY
146  if ((*scoring_previous)[st_previous] != -INFINITY){
147  viterbi_temp = getTransition((*hmm)[st_previous], st_current , position) + emission + (*scoring_previous)[st_previous];
148 
149  if (viterbi_temp > (*scoring_current)[st_current]){
150  (*scoring_current)[st_current] = viterbi_temp;
151  (*traceback_table)[position][st_current] = st_previous;
152  }
153 
154  next_states |= (*(*hmm)[st_current]->getTo());
155  }
156  }
157  }
158  }
159 
160  //TODO: Calculate ending and set the final viterbi and traceback pointer
161  //Swap current and previous viterbi scores
162  scoring_previous->assign(state_size,-INFINITY);
163  swap_ptr = scoring_previous;
164  scoring_previous = scoring_current;
165  scoring_current = swap_ptr;
166 
167 
168  //Calculate ending viterbi score and traceback from END state
169  for(size_t st_previous = 0; st_previous < state_size ;++st_previous){
170  if ((*scoring_previous)[st_previous] > -INFINITY){
171  viterbi_temp = (*scoring_previous)[st_previous] + (*hmm)[st_previous]->getEndTrans();
172 
173  if (viterbi_temp > ending_viterbi_score){
174  ending_viterbi_score = viterbi_temp;
175  ending_viterbi_tb = st_previous;
176  }
177  }
178  }
179 
180  delete scoring_previous;
181  delete scoring_current;
182  scoring_previous = NULL;
183  scoring_current = NULL;
184  }
185 
186 
187  void trellis::naive_viterbi(model* h, sequences* sqs){
188  //Initialize the table
189  hmm = h;
190  seqs = sqs;
191  seq_size = seqs->getLength();
192  state_size = hmm->state_size();
193  exDef_defined = seqs->exDefDefined();
194 
195  //TODO: determine which model and chose the type of algorithm to use;
196  naive_viterbi();
197 
198  }
199 
200 
201  void trellis::naive_viterbi(){
202  traceback_table = new(std::nothrow) int_2D(seq_size, std::vector<int16_t>(state_size,-1));
203  dbl_viterbi_score = new (std::nothrow) double_2D(seq_size, std::vector<double>(state_size,-INFINITY));
204 
205  double emission(-INFINITY);
206  double viterbi_temp(-INFINITY);
207  double trans(-INFINITY);
208  double previous(-INFINITY);
209  bool exDef_position(false);
210  ending_viterbi_tb = -1;
211  ending_viterbi_score = -INFINITY;
212 
213  state* init = hmm->getInitial();
214 
215  //Calculate from Initial states
216  for(size_t st = 0; st < state_size; ++st){
217  viterbi_temp = (*hmm)[st]->get_emission_prob(*seqs,0) + getTransition(init, st, 0);
218  (*dbl_viterbi_score)[0][st]=viterbi_temp;
219  (*traceback_table)[0][st]=-1;
220  }
221 
222  //Calculate Forward for all states
223  for (size_t position = 1 ; position < seq_size ; ++position){
224 
225  if (exDef_defined){
226  exDef_position = seqs->exDefDefined(position);
227  }
228 
229  for (size_t st_current = 0; st_current < state_size; ++st_current){
230  //Calc emissions
231  emission = (*hmm)[st_current]->get_emission_prob(*seqs, position);
232 
233  if (exDef_defined && exDef_position){
234  emission += seqs->getWeight(position, st_current);
235  }
236 
237  if (emission == -INFINITY){
238  continue;
239  }
240 
241  for (size_t st_previous = 0; st_previous < state_size; ++st_previous){
242  previous = (*dbl_viterbi_score)[position-1][st_previous];
243  if (previous == -INFINITY){
244  continue;
245  }
246 
247  trans = getTransition(hmm->getState(st_previous), st_current, position);
248 
249  if (trans !=-INFINITY){
250  viterbi_temp = emission + trans + previous;
251 
252  if (viterbi_temp > (*dbl_viterbi_score)[position][st_current]){
253  (*dbl_viterbi_score)[position][st_current] = viterbi_temp;
254  (*traceback_table)[position][st_current] = st_previous;
255  }
256  }
257  }
258  }
259  }
260 
261 
262  //Calculate Ending Transition
263  ending_viterbi_tb = -1;
264  for (size_t st_previous = 0; st_previous < state_size; ++st_previous){
265  if ((*hmm)[st_previous]->getEndTrans() != -INFINITY){
266  if ((*dbl_viterbi_score)[seq_size-1][st_previous] != -INFINITY){
267  viterbi_temp = (*dbl_viterbi_score)[seq_size-1][st_previous] + (*hmm)[st_previous]->getEndTrans();
268  if (viterbi_temp > ending_viterbi_score){
269  ending_viterbi_score = viterbi_temp;
270  ending_viterbi_tb = st_previous;
271  }
272  }
273  }
274  }
275  return;
276  }
277 
278 
279 
280 
281 
282 
283 // void trellis::viterbi(){
284 //
285 // //Initialize the traceback table
286 // if (traceback_table != NULL){
287 // delete traceback_table;
288 // }
289 //
290 // traceback_table = new int_2D(seq_size,std::vector<int16_t> (state_size,-1));
291 // scoring_previous = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
292 // scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
293 //
294 // if (scoring_previous == NULL || scoring_current == NULL || traceback_table == NULL){
295 // std::cerr << "Can't allocate Viterbi score and traceback table. OUT OF MEMORY" << std::endl;
296 // exit(2);
297 // }
298 //
299 //
300 // std::bitset<STATE_MAX> next_states;
301 // std::bitset<STATE_MAX> current_states;
302 //
303 // double viterbi_temp(-INFINITY);
304 // double emission(-INFINITY);
305 // bool exDef_position(false);
306 //
307 //
308 // //If model is not a basic model, then we need to initialize the explicit duration vector
309 // //bool extend_duration keeps track of whether the transition to same state was selected.
310 // if (!hmm->isBasic()){
311 // explicit_duration_current = new(std::nothrow) std::vector<size_t>(state_size,0);
312 // explicit_duration_previous= new(std::nothrow) std::vector<size_t>(state_size,0);
313 // }
314 // bool extend_duration(false);
315 // std::vector<bool>* duration = hmm->get_explicit();
316 //
317 // state* init = hmm->getInitial();
318 //
319 // std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
320 // std::bitset<STATE_MAX>* from_trans(NULL);
321 //
322 // //Calculate Viterbi from transitions from INIT (initial) state
323 // for(size_t i = 0; i < state_size; ++i){
324 // if ((*initial_to)[i]){ //if the bitset is set (meaning there is a transition to this state), calculate the viterbi
325 //
326 // viterbi_temp = (*hmm)[i]->get_emission_prob(*seqs,0) + getTransition(init, i, 0);
327 //
328 // if (viterbi_temp > -INFINITY){
329 // if ((*scoring_current)[i] < viterbi_temp){
330 // (*scoring_current)[i] = viterbi_temp;
331 // }
332 // next_states |= (*(*hmm)[i]->getTo());
333 // }
334 // }
335 // }
336 //
337 // // for(size_t i=0; i < state_size; ++i){
338 // // std::cout << "Position: 0" << std::endl;
339 // // std::cout << exp((*viterbi_current)[i]) << std::endl;
340 // // }
341 // //
342 //
343 // for(size_t position = 1; position < seq_size ; ++position ){
344 //
345 // //Swap current and previous viterbi scores
346 // scoring_previous->assign(state_size,-INFINITY);
347 // swap_ptr = scoring_previous;
348 // scoring_previous = scoring_current;
349 // scoring_current = swap_ptr;
350 //
351 // //Swap current_states and next states sets
352 //
353 // current_states.reset();
354 // current_states |= next_states;
355 // next_states.reset();
356 //
357 // if (exDef_defined){
358 // exDef_position = seqs->exDefDefined(position);
359 // }
360 //
361 // //TODO: Check use of external definitions below.
362 //
363 // std::cout << "\nPosition:\t" << position << "\n";
364 // // std::cout << "Letter:\t" << seqs->seqValue(0, position) << std::endl;
365 //
366 // for (size_t i = 0; i < state_size; ++i){ //i is current state that emits value
367 // if (!current_states[i]){
368 // continue;
369 // }
370 //
371 // //current_state = (*hmm)[i];
372 // //emission = current_state->get_emission(*seqs,position);
373 // emission = (*hmm)[i]->get_emission_prob(*seqs, position);
374 //
375 //
376 // // std::cout << "State Emission:\t" << i << "\t" << exp(emission) << std::endl;
377 //
378 // if (exDef_defined && exDef_position){
379 // emission += seqs->getWeight(position, i);
380 // }
381 //
382 // from_trans = (*hmm)[i]->getFrom();
383 //
384 // for (size_t j = 0; j < state_size ; ++j) { //j is previous state
385 // if (!(*from_trans)[j]){
386 // continue;
387 // }
388 //
389 // if ((*scoring_previous)[j] != -INFINITY){
390 // viterbi_temp = getTransition((*hmm)[j], i , position) + emission + (*scoring_previous)[j];
391 //
392 // std::cout << exp(getTransition((*hmm)[j],i,position)) << std::endl;
393 //
394 // // std::cout << "Temp Viterbi:\tTransFrom: "<< j << "\tto\t" << i << "\t" << viterbi_temp / log(2) << std::endl;
395 //
396 //
397 // if (viterbi_temp > (*scoring_current)[i]){
398 // //If transition is from same to same then if it is
399 // //explicit duration we'll need to change
400 // extend_duration = (i==j) ? true : false;
401 //
402 // (*scoring_current)[i] = viterbi_temp;
403 // (*traceback_table)[position][i] = j;
404 // }
405 //
406 // next_states |= (*(*hmm)[i]->getTo());
407 // }
408 // }
409 //
410 // //If explicit durations vector defined, and transition from-to same
411 // //then we'll increment the value. Otherwise, set to zero;
412 // if (explicit_duration_current){
413 // if (extend_duration && (*duration)[i]){
414 // (*explicit_duration_current)[i]=(*explicit_duration_previous)[i]+1;
415 // extend_duration=false;
416 // }
417 // else{
418 // (*explicit_duration_current)[i]=0;
419 // }
420 // }
421 // }
422 //
423 // if(explicit_duration_current){
424 // swap_ptr_duration = explicit_duration_previous;
425 // explicit_duration_previous = explicit_duration_current;
426 // explicit_duration_current= swap_ptr_duration;
427 // explicit_duration_current->assign(state_size,0);
428 // }
429 //
430 // }
431 //
432 // //TODO: Calculate ending and set the final viterbi and traceback pointer
433 // //Swap current and previous viterbi scores
434 // scoring_previous->assign(state_size,-INFINITY);
435 // swap_ptr = scoring_previous;
436 // scoring_previous = scoring_current;
437 // scoring_current = swap_ptr;
438 //
439 // for(size_t i = 0; i < state_size ;++i){
440 // if ((*scoring_previous)[i] > -INFINITY){
441 // viterbi_temp = (*scoring_previous)[i] + (*hmm)[i]->getEndTrans();
442 //
443 // if (viterbi_temp > ending_viterbi_score){
444 // ending_viterbi_score = viterbi_temp;
445 // ending_viterbi_tb = i;
446 // }
447 // }
448 // }
449 //
450 // delete scoring_previous;
451 // delete scoring_current;
452 // scoring_previous = NULL;
453 // scoring_current = NULL;
454 // }
455 
456 
457  //TODO: Need to test and finalize these algorithms perform
458  /* Need to simplify the basic calls so that it checks model and chooses the
459  algorithm to perform
460 
461  Need to establish duration algorithms for forward/backward that first to viterbi
462  to calculate the transition probability.
463 
464  */
465 
466 
467  //! Sparse Complex Viterbi
468  //! Stores the transition duration probabilities in a hashmap (memory efficient, slower)
469  //! The duratio probabilities can then be used in forward and backward algorithms.
470  void trellis::sparse_complex_viterbi(){
471 
472  //Initialize the traceback table
473  if (traceback_table != NULL){
474  delete traceback_table;
475  }
476 
477  traceback_table = new int_2D(seq_size,std::vector<int16_t> (state_size,1));
478  scoring_previous = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
479  scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
480 
481  if (scoring_previous == NULL || scoring_current == NULL || traceback_table == NULL){
482  std::cerr << "Can't allocate Viterbi score and traceback table. OUT OF MEMORY" << std::endl;
483  exit(2);
484  }
485 
486 
487  std::bitset<STATE_MAX> next_states;
488  std::bitset<STATE_MAX> current_states;
489 
490  double viterbi_temp(-INFINITY);
491  double emission(-INFINITY);
492  bool exDef_position(false);
493  double transition_prob(-INFINITY);
494  ending_viterbi_score = -INFINITY;
495  ending_viterbi_tb = -1;
496 
497 
498  //If model is not a basic model, then we need to initialize the explicit duration vector
499  //bool extend_duration keeps track of whether the transition to same state was selected.
500  if (!hmm->isBasic()){
501  explicit_duration_current = new(std::nothrow) std::vector<size_t>(state_size,0);
502  explicit_duration_previous= new(std::nothrow) std::vector<size_t>(state_size,0);
503  }
504 
505  bool extend_duration(false);
506 
507  // Get list of States with explicit duration
508  std::vector<bool>* duration = hmm->get_explicit();
509 
510  //Initialize Duration storage table
511  complex_transitions = new std::vector<std::vector< std::map<uint16_t,double>* >* > (state_size,NULL);
512  for(size_t i=0; i < state_size; i++){
513  if((*duration)[i]){
514  (*complex_transitions)[i] = new std::vector<std::map<uint16_t,double>* > (seq_size, NULL);
515  }
516  }
517 
518  state* init = hmm->getInitial();
519 
520  std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
521  std::bitset<STATE_MAX>* from_trans(NULL);
522 
523  //Calculate Viterbi from transitions from INIT (initial) state
524  for(size_t i = 0; i < state_size; ++i){
525  if ((*initial_to)[i]){ //if the bitset is set (meaning there is a transition to this state), calculate the viterbi
526 
527  //Transitions here are guarenteed to be standard from the initial state
528  viterbi_temp = (*hmm)[i]->get_emission_prob(*seqs,0) + getTransition(init, i, 0);
529 
530  if (viterbi_temp > -INFINITY){
531  if ((*scoring_current)[i] < viterbi_temp){
532  (*scoring_current)[i] = viterbi_temp;
533  }
534  next_states |= (*(*hmm)[i]->getTo());
535  }
536  }
537  }
538 
539 
540  for(size_t position = 1; position < seq_size ; ++position ){
541 
542  //Swap current and previous viterbi scores
543  scoring_previous->assign(state_size,-INFINITY);
544  swap_ptr = scoring_previous;
545  scoring_previous = scoring_current;
546  scoring_current = swap_ptr;
547 
548  //Swap current_states and next states sets
549 
550  current_states.reset();
551  current_states |= next_states;
552  next_states.reset();
553 
554  if (exDef_defined){
555  exDef_position = seqs->exDefDefined(position);
556  }
557 
558  //TODO: Check use of external definitions below.
559 
560  //std::cout << "\nPosition:\t" << position << "\n";
561  // std::cout << "Letter:\t" << seqs->seqValue(0, position) << std::endl;
562 
563  for (size_t st_current = 0; st_current < state_size; ++st_current){ //i is current state that emits value
564  if (!current_states[st_current]){
565  continue;
566  }
567 
568  emission = (*hmm)[st_current]->get_emission_prob(*seqs, position);
569 
570 
571  //Check External definitions
572  if (exDef_defined && exDef_position){
573  emission += seqs->getWeight(position, st_current);
574  }
575 
576  //Get list of state that transition to current state
577  from_trans = (*hmm)[st_current]->getFrom();
578 
579  for (size_t st_previous = 0; st_previous < state_size ; ++st_previous) { //j is previous state
580  if (!(*from_trans)[st_previous]){ //if transition not possible next
581  continue;
582  }
583 
584  if ((*scoring_previous)[st_previous] != -INFINITY){
585 
586  transition_prob = getTransition((*hmm)[st_previous], st_current , position);
587 
588  if ((*duration)[st_previous]){
589  if ((*(*complex_transitions)[st_current])[position] == NULL){
590  (*(*complex_transitions)[st_current])[position] = new std::map<uint16_t,double>;
591  }
592 
593  (*(*(*complex_transitions)[st_current])[position])[st_previous] = transition_prob;
594  }
595 
596  viterbi_temp = transition_prob + emission + (*scoring_previous)[st_previous];
597 
598 
599  if (viterbi_temp > (*scoring_current)[st_current]){
600  //If transition is from same to same then if it is
601  //explicit duration we'll need to change
602  extend_duration = (st_current==st_previous) ? true : false;
603 
604  (*scoring_current)[st_current] = viterbi_temp;
605  (*traceback_table)[position][st_current] = st_previous;
606  }
607 
608  next_states |= (*(*hmm)[st_current]->getTo());
609  }
610  }
611 
612  //If explicit durations vector defined, and transition from-to same
613  //then we'll increment the value. Otherwise, set to zero;
614  if (explicit_duration_current){
615  if (extend_duration && (*duration)[st_current]){
616  if ((*explicit_duration_previous)[st_current]==0){
617  (*explicit_duration_current)[st_current]=2;
618  }
619  else{
620  (*explicit_duration_current)[st_current]=(*explicit_duration_previous)[st_current]+1;
621  }
622  extend_duration=false;
623  }
624  else{
625  (*explicit_duration_current)[st_current]=0;
626  }
627  }
628  }
629 
630  if(explicit_duration_current){
631  swap_ptr_duration = explicit_duration_previous;
632  explicit_duration_previous = explicit_duration_current;
633  explicit_duration_current= swap_ptr_duration;
634  explicit_duration_current->assign(state_size,0);
635  }
636 
637  }
638 
639  //TODO: Calculate ending and set the final viterbi and traceback pointer
640  //Swap current and previous viterbi scores
641  scoring_previous->assign(state_size,-INFINITY);
642  swap_ptr = scoring_previous;
643  scoring_previous = scoring_current;
644  scoring_current = swap_ptr;
645 
646  for(size_t st_previous = 0; st_previous < state_size ;++st_previous){
647  if ((*scoring_previous)[st_previous] > -INFINITY){
648  viterbi_temp = (*scoring_previous)[st_previous] + (*hmm)[st_previous]->getEndTrans();
649 
650  if (viterbi_temp > ending_viterbi_score){
651  ending_viterbi_score = viterbi_temp;
652  ending_viterbi_tb = st_previous;
653  }
654  }
655  }
656 
657  delete scoring_previous;
658  delete scoring_current;
659  scoring_previous = NULL;
660  scoring_current = NULL;
661  }
662 
663  //!Store transition in a table (more memory, but faster than sparse complex)
664  //!Need testing and more develepment;
665  void trellis::fast_complex_viterbi(){
666 
667  //Initialize the traceback table
668  if (traceback_table != NULL){
669  delete traceback_table;
670  }
671 
672  traceback_table = new int_2D(seq_size,std::vector<int16_t> (state_size,-1));
673  scoring_previous = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
674  scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
675 
676  if (scoring_previous == NULL || scoring_current == NULL || traceback_table == NULL){
677  std::cerr << "Can't allocate Viterbi score and traceback table. OUT OF MEMORY" << std::endl;
678  exit(2);
679  }
680 
681 
682  std::bitset<STATE_MAX> next_states;
683  std::bitset<STATE_MAX> current_states;
684 
685  double viterbi_temp(-INFINITY);
686  double emission(-INFINITY);
687  bool exDef_position(false);
688  ending_viterbi_tb = -1;
689  ending_viterbi_score = -INFINITY;
690 
691 
692  //If model is not a basic model, then we need to initialize the explicit duration vector
693  //bool extend_duration keeps track of whether the transition to same state was selected.
694  if (!hmm->isBasic()){
695  explicit_duration_current = new(std::nothrow) std::vector<size_t>(state_size,0);
696  explicit_duration_previous= new(std::nothrow) std::vector<size_t>(state_size,0);
697  }
698  bool extend_duration(false);
699  std::vector<bool>* duration = hmm->get_explicit();
700 
701  state* init = hmm->getInitial();
702 
703  std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
704  std::bitset<STATE_MAX>* from_trans(NULL);
705 
706  //Calculate Viterbi from transitions from INIT (initial) state
707  for(size_t i = 0; i < state_size; ++i){
708  if ((*initial_to)[i]){ //if the bitset is set (meaning there is a transition to this state), calculate the viterbi
709 
710  viterbi_temp = (*hmm)[i]->get_emission_prob(*seqs,0) + getTransition(init, i, 0);
711 
712  if (viterbi_temp > -INFINITY){
713  if ((*scoring_current)[i] < viterbi_temp){
714  (*scoring_current)[i] = viterbi_temp;
715  }
716  next_states |= (*(*hmm)[i]->getTo());
717 // if ((*duration)[i]){
718 // (*explicit_duration_previous)[i]=1;
719 // }
720  }
721  }
722  }
723 
724 
725  for(size_t position = 1; position < seq_size ; ++position ){
726 
727  //Swap current and previous viterbi scores
728  scoring_previous->assign(state_size,-INFINITY);
729  swap_ptr = scoring_previous;
730  scoring_previous = scoring_current;
731  scoring_current = swap_ptr;
732 
733  //Swap current_states and next states sets
734 
735  current_states.reset();
736  current_states |= next_states;
737  next_states.reset();
738 
739  if (exDef_defined){
740  exDef_position = seqs->exDefDefined(position);
741  }
742 
743  //TODO: Check use of external definitions below.
744 
745 // std::cout << "\nPosition:\t" << position << "\n";
746 // std::cout << "Letter:\t" << seqs->seqValue(0, position)+1 << std::endl;
747 
748  for (size_t i = 0; i < state_size; ++i){ //i is current state that emits value
749  if (!current_states[i]){
750  continue;
751  }
752 
753  //current_state = (*hmm)[i];
754  //emission = current_state->get_emission(*seqs,position);
755  emission = (*hmm)[i]->get_emission_prob(*seqs, position);
756 
757 
758 // std::cout << "State Emission:\t" << i << "\t" << exp(emission) << std::endl;
759 
760  if (exDef_defined && exDef_position){
761  emission += seqs->getWeight(position, i);
762  }
763 
764  from_trans = (*hmm)[i]->getFrom();
765 
766  for (size_t j = 0; j < state_size ; ++j) { //j is previous state
767  if (!(*from_trans)[j]){
768  continue;
769  }
770 
771  if ((*scoring_previous)[j] != -INFINITY){
772  viterbi_temp = getTransition((*hmm)[j], i , position) + emission + (*scoring_previous)[j];
773 
774 // std::cout << "TransFrom: "<< j << "\tto\t" << i << "\t" << exp(getTransition((*hmm)[j],i,position)) << std::endl;
775 //
776 // std::cout << "Temp Viterbi:\t" << viterbi_temp << std::endl;
777 
778 
779  if (viterbi_temp > (*scoring_current)[i]){
780  //If transition is from same to same then if it is
781  //explicit duration we'll need to change
782  extend_duration = (i==j) ? true : false;
783 
784  (*scoring_current)[i] = viterbi_temp;
785  (*traceback_table)[position][i] = j;
786  }
787 
788  next_states |= (*(*hmm)[i]->getTo());
789  }
790  }
791 
792  //If explicit durations vector defined, and transition from-to same
793  //then we'll increment the value. Otherwise, set to zero;
794  if (explicit_duration_current){
795  if (extend_duration && (*duration)[i]){
796  if ((*explicit_duration_previous)[i]==0){
797  (*explicit_duration_current)[i]=2;
798  }
799  else{
800  (*explicit_duration_current)[i]=(*explicit_duration_previous)[i]+1;
801  }
802  extend_duration=false;
803  }
804  else{
805  (*explicit_duration_current)[i]=0;
806  }
807  }
808  }
809 
810  if(explicit_duration_current){
811  swap_ptr_duration = explicit_duration_previous;
812  explicit_duration_previous = explicit_duration_current;
813  explicit_duration_current= swap_ptr_duration;
814  explicit_duration_current->assign(state_size,0);
815  }
816 
817  }
818 
819  //TODO: Calculate ending and set the final viterbi and traceback pointer
820  //Swap current and previous viterbi scores
821  scoring_previous->assign(state_size,-INFINITY);
822  swap_ptr = scoring_previous;
823  scoring_previous = scoring_current;
824  scoring_current = swap_ptr;
825 
826  for(size_t i = 0; i < state_size ;++i){
827  if ((*scoring_previous)[i] > -INFINITY){
828  viterbi_temp = (*scoring_previous)[i] + (*hmm)[i]->getEndTrans();
829 
830  if (viterbi_temp > ending_viterbi_score){
831  ending_viterbi_score = viterbi_temp;
832  ending_viterbi_tb = i;
833  }
834  }
835  }
836 
837  delete scoring_previous;
838  delete scoring_current;
839  scoring_previous = NULL;
840  scoring_current = NULL;
841  }
842 
843 
844 }
845 
Generated on Tue Jul 30 2013 13:23:12 for StochHMM by   doxygen 1.8.1