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_forward.cpp
Go to the documentation of this file.
1 //
2 // stoch_forward.cpp
3 // StochHMM
4 //
5 // Created by Paul Lott on 2/6/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_forward(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  stochastic_forward();
21  }
22 
23 
24  void trellis::stochastic_forward(){
25  //if (hmm->isBasic()){
26  simple_stochastic_forward();
27  //}
28  }
29 
30  void trellis::simple_stochastic_forward(){
31 
32 // if (!hmm->isBasic()){
33 // std::cerr << "Model isn't a simple/basic HMM. Use complex algorithms\n";
34 // return;
35 // }
36 
37 
38  scoring_current = new (std::nothrow) std::vector<double> (state_size,-INFINITY);
39  scoring_previous= new (std::nothrow) std::vector<double> (state_size,-INFINITY);
40  stochastic_table = new (std::nothrow) stochTable(seq_size);
41 
42  if (scoring_current == NULL || scoring_previous == NULL || stochastic_table == NULL){
43  std::cerr << "Can't allocate forward score table. OUT OF MEMORY" << std::endl;
44  exit(2);
45  }
46 
47  std::bitset<STATE_MAX> next_states;
48  std::bitset<STATE_MAX> current_states;
49 
50  double forward_temp(-INFINITY);
51  double emission(-INFINITY);
52  bool exDef_position(false);
53 
54  state* init = hmm->getInitial();
55 
56  std::bitset<STATE_MAX>* initial_to = hmm->getInitialTo();
57  std::bitset<STATE_MAX>* from_trans(NULL);
58 
59 
60  // std::cout << "Position: 0" << std::endl;
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  forward_temp = (*hmm)[st]->get_emission_prob(*seqs,0) + getTransition(init, st, 0);
66 
67  if (forward_temp > -INFINITY){
68  (*scoring_current)[st] = forward_temp;
69  next_states |= (*(*hmm)[st]->getTo());
70  }
71  }
72  }
73 
74 
75  for(size_t position = 1; position < seq_size ; ++position ){
76 
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 
85  //Swap current_states and next states sets
86  current_states.reset();
87  current_states |= next_states;
88  next_states.reset();
89 
90  if (exDef_defined){
91  exDef_position = seqs->exDefDefined(position);
92  }
93 
94  // std::cout << "\nPosition: " << position << std::endl;
95 
96  for (size_t st_current = 0; st_current < state_size; ++st_current){ //i is current state that emits value
97  if (!current_states[st_current]){
98  continue;
99  }
100 
101  emission = (*hmm)[st_current]->get_emission_prob(*seqs, position);
102 
103  if (exDef_defined && exDef_position){
104  emission += seqs->getWeight(position, st_current);
105  }
106 
107  from_trans = (*hmm)[st_current]->getFrom();
108 
109  for (size_t st_previous = 0; st_previous < state_size ; ++st_previous) { //j is previous state
110  if (!(*from_trans)[st_previous]){
111  continue;
112  }
113 
114  if ((*scoring_previous)[st_previous] != -INFINITY){
115  forward_temp = (*scoring_previous)[st_previous] + emission + getTransition((*hmm)[st_previous], st_current , position);
116 
117  if (forward_temp == -INFINITY){
118  continue;
119  }
120 
121  stochastic_table->push(position-1,st_current,st_previous,forward_temp);
122 
123  if ((*scoring_current)[st_current] == -INFINITY){
124  (*scoring_current)[st_current] = forward_temp;
125  }
126  else{
127  (*scoring_current)[st_current] = addLog(forward_temp, (*scoring_current)[st_current]);
128  }
129 
130  next_states |= (*(*hmm)[st_current]->getTo());
131  }
132  }
133  }
134  }
135 
136  //Swap current and previous scores
137  scoring_previous->assign(state_size,-INFINITY);
138  swap_ptr = scoring_previous;
139  scoring_previous = scoring_current;
140  scoring_current = swap_ptr;
141 
142  ending_forward_prob = -INFINITY;
143  for(size_t st_previous = 0; st_previous < state_size ;++st_previous){
144  if ((*scoring_previous)[st_previous] != -INFINITY){
145  forward_temp = (*scoring_previous)[st_previous] + (*hmm)[st_previous]->getEndTrans();
146 
147  if (forward_temp == -INFINITY){
148  continue;
149  }
150  stochastic_table->push(seq_size-1,SIZE_MAX, st_previous,forward_temp);
151 
152  if (forward_temp > -INFINITY){
153  if (ending_forward_prob == -INFINITY){
154  ending_forward_prob = forward_temp;
155  }
156  else{
157  ending_forward_prob = addLog(ending_forward_prob,forward_temp);
158  }
159  }
160  }
161  }
162 
163  stochastic_table->finalize();
164  //stochastic_table->print();
165 
166  delete scoring_previous;
167  delete scoring_current;
168  scoring_previous = NULL;
169  scoring_current = NULL;
170 
171  }
172 
173 
174  void trellis::naive_stochastic_forward(){
175  dbl_forward_score = new (std::nothrow) double_2D(seq_size, std::vector<double>(state_size,-INFINITY));
176  stochastic_table = new (std::nothrow) stochTable(seq_size);
177 
178 
179  if (dbl_forward_score == NULL){
180  std::cerr << "Can't allocate Forward table and traceback table. OUT OF MEMORY\t" << __FUNCTION__ << std::endl;
181  exit(2);
182  }
183 
184  double emission(-INFINITY);
185  double forward_temp(-INFINITY);
186  double trans(-INFINITY);
187  double previous(-INFINITY);
188  bool exDef_position(false);
189 
190  state* init = hmm->getInitial();
191 
192  //Calculate from Initial states
193  for(size_t st = 0; st < state_size; ++st){
194  forward_temp = (*hmm)[st]->get_emission_prob(*seqs,0) + getTransition(init, st, 0);
195  (*dbl_forward_score)[0][st]=forward_temp;
196  }
197 
198  //Calculate Forward for all states
199  for (size_t position = 1 ; position < seq_size ; ++position){
200  for (size_t st_current = 0; st_current < state_size; ++st_current){
201  //Calc emissions
202  emission = (*hmm)[st_current]->get_emission_prob(*seqs, position);
203 
204  if (exDef_defined && exDef_position){
205  emission += seqs->getWeight(position, st_current);
206  }
207 
208  if (emission == -INFINITY){
209  continue;
210  }
211 
212  for (size_t st_previous = 0; st_previous < state_size; ++st_previous){
213  previous = (*dbl_forward_score)[position-1][st_previous];
214 
215  if (previous == -INFINITY){
216  continue;
217  }
218 
219  trans = getTransition(hmm->getState(st_previous), st_current, position);
220 
221 
222  if (trans !=-INFINITY ){
223  forward_temp = previous + emission + trans;
224 
225  //Save partial value to stochastic table
226  stochastic_table->push(position-1,st_current,st_previous,forward_temp);
227 
228  if ((*dbl_forward_score)[position][st_current]==-INFINITY){
229  (*dbl_forward_score)[position][st_current] = forward_temp;
230  }
231  else{
232  (*dbl_forward_score)[position][st_current] = addLog(forward_temp, (*dbl_forward_score)[position][st_current]);
233  }
234  }
235  }
236  }
237  }
238 
239 
240  //Calculate Ending Transition
241  ending_forward_prob = -INFINITY;
242  for (size_t st_previous = 0; st_previous < state_size; ++st_previous){
243 
244  if ((*hmm)[st_previous]->getEndTrans() != -INFINITY){
245 
246  if ((*dbl_forward_score)[seq_size-1][st_previous] != -INFINITY){
247  forward_temp = (*dbl_forward_score)[seq_size-1][st_previous] + (*hmm)[st_previous]->getEndTrans();
248  stochastic_table->push(seq_size-1,SIZE_MAX, st_previous,forward_temp);
249 
250  ending_forward_prob = addLog(forward_temp, ending_forward_prob);
251  }
252  }
253  }
254 
255  stochastic_table->finalize();
256  //stochastic_table->print();
257 
258  return;
259  }
260 
261 }
Generated on Tue Jul 30 2013 13:23:12 for StochHMM by   doxygen 1.8.1