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