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