SHOGUN  v3.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HMM.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 1999-2009 Soeren Sonnenburg
8  * Written (W) 1999-2008 Gunnar Raetsch
9  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
10  */
11 #include <shogun/distributions/HMM.h>
12 #include <shogun/mathematics/Math.h>
13 #include <shogun/io/SGIO.h>
14 #include <shogun/lib/config.h>
15 #include <shogun/lib/Signal.h>
16 #include <shogun/base/Parallel.h>
17 #include <shogun/features/StringFeatures.h>
18 #include <shogun/features/Alphabet.h>
19 
20 #include <stdlib.h>
21 #include <stdio.h>
22 #include <time.h>
23 #include <ctype.h>
24 
25 #define VAL_MACRO log((default_value == 0) ? (CMath::random(MIN_RAND, MAX_RAND)) : default_value)
26 #define ARRAY_SIZE 65336
27 
28 using namespace shogun;
29 
31 // Construction/Destruction
33 
34 const int32_t CHMM::GOTN= (1<<1);
35 const int32_t CHMM::GOTM= (1<<2);
36 const int32_t CHMM::GOTO= (1<<3);
37 const int32_t CHMM::GOTa= (1<<4);
38 const int32_t CHMM::GOTb= (1<<5);
39 const int32_t CHMM::GOTp= (1<<6);
40 const int32_t CHMM::GOTq= (1<<7);
41 
42 const int32_t CHMM::GOTlearn_a= (1<<1);
43 const int32_t CHMM::GOTlearn_b= (1<<2);
44 const int32_t CHMM::GOTlearn_p= (1<<3);
45 const int32_t CHMM::GOTlearn_q= (1<<4);
46 const int32_t CHMM::GOTconst_a= (1<<5);
47 const int32_t CHMM::GOTconst_b= (1<<6);
48 const int32_t CHMM::GOTconst_p= (1<<7);
49 const int32_t CHMM::GOTconst_q= (1<<8);
50 
51 enum E_STATE
52 {
53  INITIAL,
54  ARRAYs,
55  GET_N,
56  GET_M,
57  GET_a,
58  GET_b,
59  GET_p,
60  GET_q,
61  GET_learn_a,
62  GET_learn_b,
63  GET_learn_p,
64  GET_learn_q,
65  GET_const_a,
66  GET_const_b,
67  GET_const_p,
68  GET_const_q,
69  COMMENT,
70  END
71 };
72 
73 
74 #ifdef FIX_POS
75 const char Model::FIX_DISALLOWED=0 ;
76 const char Model::FIX_ALLOWED=1 ;
77 const char Model::FIX_DEFAULT=-1 ;
78 const float64_t Model::DISALLOWED_PENALTY=CMath::ALMOST_NEG_INFTY ;
79 #endif
80 
81 Model::Model()
82 {
83  const_a=SG_MALLOC(int, ARRAY_SIZE);
84  const_b=SG_MALLOC(int, ARRAY_SIZE);
85  const_p=SG_MALLOC(int, ARRAY_SIZE);
86  const_q=SG_MALLOC(int, ARRAY_SIZE);
87  const_a_val=SG_MALLOC(float64_t, ARRAY_SIZE);
88  const_b_val=SG_MALLOC(float64_t, ARRAY_SIZE);
89  const_p_val=SG_MALLOC(float64_t, ARRAY_SIZE);
90  const_q_val=SG_MALLOC(float64_t, ARRAY_SIZE);
91 
92 
93  learn_a=SG_MALLOC(int, ARRAY_SIZE);
94  learn_b=SG_MALLOC(int, ARRAY_SIZE);
95  learn_p=SG_MALLOC(int, ARRAY_SIZE);
96  learn_q=SG_MALLOC(int, ARRAY_SIZE);
97 
98 #ifdef FIX_POS
99  fix_pos_state = SG_MALLOC(char, ARRAY_SIZE);
100 #endif
101  for (int32_t i=0; i<ARRAY_SIZE; i++)
102  {
103  const_a[i]=-1 ;
104  const_b[i]=-1 ;
105  const_p[i]=-1 ;
106  const_q[i]=-1 ;
107  const_a_val[i]=1.0 ;
108  const_b_val[i]=1.0 ;
109  const_p_val[i]=1.0 ;
110  const_q_val[i]=1.0 ;
111  learn_a[i]=-1 ;
112  learn_b[i]=-1 ;
113  learn_p[i]=-1 ;
114  learn_q[i]=-1 ;
115 #ifdef FIX_POS
116  fix_pos_state[i] = FIX_DEFAULT ;
117 #endif
118  } ;
119 }
120 
121 Model::~Model()
122 {
123  SG_FREE(const_a);
124  SG_FREE(const_b);
125  SG_FREE(const_p);
126  SG_FREE(const_q);
127  SG_FREE(const_a_val);
128  SG_FREE(const_b_val);
129  SG_FREE(const_p_val);
130  SG_FREE(const_q_val);
131 
132  SG_FREE(learn_a);
133  SG_FREE(learn_b);
134  SG_FREE(learn_p);
135  SG_FREE(learn_q);
136 
137 #ifdef FIX_POS
138  SG_FREE(fix_pos_state);
139 #endif
140 
141 }
142 
143 CHMM::CHMM()
144 {
145  N=0;
146  M=0;
147  model=NULL;
148  status=false;
149  p_observations=NULL;
150  trans_list_forward_cnt=NULL;
151  trans_list_backward_cnt=NULL;
152  trans_list_forward=NULL;
153  trans_list_backward=NULL;
154  trans_list_forward_val=NULL;
155  iterations=150;
156  epsilon=1e-4;
157  conv_it=5;
158  path=NULL;
159  arrayN1=NULL;
160  arrayN2=NULL;
161  reused_caches=false;
162  transition_matrix_a=NULL;
163  observation_matrix_b=NULL;
164  initial_state_distribution_p=NULL;
165  end_state_distribution_q=NULL;
166 #ifdef USE_HMMPARALLEL_STRUCTURES
167  this->alpha_cache=NULL;
168  this->beta_cache=NULL;
169 #else
170  this->alpha_cache.table=NULL;
171  this->beta_cache.table=NULL;
172  this->alpha_cache.dimension=0;
173  this->beta_cache.dimension=0;
174 #endif
175  states_per_observation_psi=NULL;
176 }
177 
178 CHMM::CHMM(CHMM* h)
179 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
180 {
181 #ifdef USE_HMMPARALLEL_STRUCTURES
182  SG_INFO("hmm is using %i separate tables\n", parallel->get_num_threads())
183 #endif
184 
185  this->N=h->get_N();
186  this->M=h->get_M();
187  status=initialize(NULL, h->get_pseudo());
188  this->copy_model(h);
189  set_observations(h->p_observations);
190 }
191 
192 CHMM::CHMM(int32_t p_N, int32_t p_M, Model* p_model, float64_t p_PSEUDO)
193 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
194 {
195  this->N=p_N;
196  this->M=p_M;
197  model=NULL ;
198 
199 #ifdef USE_HMMPARALLEL_STRUCTURES
200  SG_INFO("hmm is using %i separate tables\n", parallel->get_num_threads())
201 #endif
202 
203  status=initialize(p_model, p_PSEUDO);
204 }
205 
206 CHMM::CHMM(
207  CStringFeatures<uint16_t>* obs, int32_t p_N, int32_t p_M,
208  float64_t p_PSEUDO)
209 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
210 {
211  this->N=p_N;
212  this->M=p_M;
213  model=NULL ;
214 
215 #ifdef USE_HMMPARALLEL_STRUCTURES
216  SG_INFO("hmm is using %i separate tables\n", parallel->get_num_threads())
217 #endif
218 
219  initialize(model, p_PSEUDO);
220  set_observations(obs);
221 }
222 
223 CHMM::CHMM(int32_t p_N, float64_t* p, float64_t* q, float64_t* a)
224 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
225 {
226  this->N=p_N;
227  this->M=0;
228  model=NULL ;
229 
230  trans_list_forward = NULL ;
231  trans_list_forward_cnt = NULL ;
232  trans_list_forward_val = NULL ;
233  trans_list_backward = NULL ;
234  trans_list_backward_cnt = NULL ;
235  trans_list_len = 0 ;
236  mem_initialized = false ;
237 
238  this->transition_matrix_a=NULL;
239  this->observation_matrix_b=NULL;
240  this->initial_state_distribution_p=NULL;
241  this->end_state_distribution_q=NULL;
242  this->p_observations=NULL;
243  this->reused_caches=false;
244 
245 #ifdef USE_HMMPARALLEL_STRUCTURES
246  this->alpha_cache=NULL;
247  this->beta_cache=NULL;
248 #else
249  this->alpha_cache.table=NULL;
250  this->beta_cache.table=NULL;
251  this->alpha_cache.dimension=0;
252  this->beta_cache.dimension=0;
253 #endif
254 
255  this->states_per_observation_psi=NULL ;
256  this->path=NULL;
257  arrayN1=NULL ;
258  arrayN2=NULL ;
259 
260  this->loglikelihood=false;
261  mem_initialized = true ;
262 
263  transition_matrix_a=a ;
264  observation_matrix_b=NULL ;
265  initial_state_distribution_p=p ;
266  end_state_distribution_q=q ;
267  transition_matrix_A=NULL ;
268  observation_matrix_B=NULL ;
269 
270 // this->invalidate_model();
271 }
272 
273 CHMM::CHMM(
274  int32_t p_N, float64_t* p, float64_t* q, int32_t num_trans,
275  float64_t* a_trans)
276 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
277 {
278  model=NULL ;
279 
280  this->N=p_N;
281  this->M=0;
282 
283  trans_list_forward = NULL ;
284  trans_list_forward_cnt = NULL ;
285  trans_list_forward_val = NULL ;
286  trans_list_backward = NULL ;
287  trans_list_backward_cnt = NULL ;
288  trans_list_len = 0 ;
289  mem_initialized = false ;
290 
291  this->transition_matrix_a=NULL;
292  this->observation_matrix_b=NULL;
293  this->initial_state_distribution_p=NULL;
294  this->end_state_distribution_q=NULL;
295  this->p_observations=NULL;
296  this->reused_caches=false;
297 
298 #ifdef USE_HMMPARALLEL_STRUCTURES
299  this->alpha_cache=NULL;
300  this->beta_cache=NULL;
301 #else
302  this->alpha_cache.table=NULL;
303  this->beta_cache.table=NULL;
304  this->alpha_cache.dimension=0;
305  this->beta_cache.dimension=0;
306 #endif
307 
308  this->states_per_observation_psi=NULL ;
309  this->path=NULL;
310  arrayN1=NULL ;
311  arrayN2=NULL ;
312 
313  this->loglikelihood=false;
314  mem_initialized = true ;
315 
316  trans_list_forward_cnt=NULL ;
317  trans_list_len = N ;
318  trans_list_forward = SG_MALLOC(T_STATES*, N);
319  trans_list_forward_val = SG_MALLOC(float64_t*, N);
320  trans_list_forward_cnt = SG_MALLOC(T_STATES, N);
321 
322  int32_t start_idx=0;
323  for (int32_t j=0; j<N; j++)
324  {
325  int32_t old_start_idx=start_idx;
326 
327  while (start_idx<num_trans && a_trans[start_idx+num_trans]==j)
328  {
329  start_idx++;
330 
331  if (start_idx>1 && start_idx<num_trans)
332  ASSERT(a_trans[start_idx+num_trans-1]<=
333  a_trans[start_idx+num_trans]);
334  }
335 
336  if (start_idx>1 && start_idx<num_trans)
337  ASSERT(a_trans[start_idx+num_trans-1]<=
338  a_trans[start_idx+num_trans]);
339 
340  int32_t len=start_idx-old_start_idx;
341  ASSERT(len>=0)
342 
343  trans_list_forward_cnt[j] = 0 ;
344 
345  if (len>0)
346  {
347  trans_list_forward[j] = SG_MALLOC(T_STATES, len);
348  trans_list_forward_val[j] = SG_MALLOC(float64_t, len);
349  }
350  else
351  {
352  trans_list_forward[j] = NULL;
353  trans_list_forward_val[j] = NULL;
354  }
355  }
356 
357  for (int32_t i=0; i<num_trans; i++)
358  {
359  int32_t from = (int32_t)a_trans[i+num_trans] ;
360  int32_t to = (int32_t)a_trans[i] ;
361  float64_t val = a_trans[i+num_trans*2] ;
362 
363  ASSERT(from>=0 && from<N)
364  ASSERT(to>=0 && to<N)
365 
366  trans_list_forward[from][trans_list_forward_cnt[from]]=to ;
367  trans_list_forward_val[from][trans_list_forward_cnt[from]]=val ;
368  trans_list_forward_cnt[from]++ ;
369  //ASSERT(trans_list_forward_cnt[from]<3000)
370  } ;
371 
372  transition_matrix_a=NULL ;
373  observation_matrix_b=NULL ;
374  initial_state_distribution_p=p ;
375  end_state_distribution_q=q ;
376  transition_matrix_A=NULL ;
377  observation_matrix_B=NULL ;
378 
379 // this->invalidate_model();
380 }
381 
382 
383 CHMM::CHMM(FILE* model_file, float64_t p_PSEUDO)
384 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
385 {
386 #ifdef USE_HMMPARALLEL_STRUCTURES
387  SG_INFO("hmm is using %i separate tables\n", parallel->get_num_threads())
388 #endif
389 
390  status=initialize(NULL, p_PSEUDO, model_file);
391 }
392 
393 CHMM::~CHMM()
394 {
395  SG_UNREF(p_observations);
396 
397  if (trans_list_forward_cnt)
398  SG_FREE(trans_list_forward_cnt);
399  if (trans_list_backward_cnt)
400  SG_FREE(trans_list_backward_cnt);
401  if (trans_list_forward)
402  {
403  for (int32_t i=0; i<trans_list_len; i++)
404  if (trans_list_forward[i])
405  SG_FREE(trans_list_forward[i]);
406  SG_FREE(trans_list_forward);
407  }
408  if (trans_list_forward_val)
409  {
410  for (int32_t i=0; i<trans_list_len; i++)
411  if (trans_list_forward_val[i])
412  SG_FREE(trans_list_forward_val[i]);
413  SG_FREE(trans_list_forward_val);
414  }
415  if (trans_list_backward)
416  {
417  for (int32_t i=0; i<trans_list_len; i++)
418  if (trans_list_backward[i])
419  SG_FREE(trans_list_backward[i]);
420  SG_FREE(trans_list_backward);
421  } ;
422 
423  free_state_dependend_arrays();
424 
425  if (!reused_caches)
426  {
427 #ifdef USE_HMMPARALLEL_STRUCTURES
428  for (int32_t i=0; i<parallel->get_num_threads(); i++)
429  {
430  SG_FREE(alpha_cache[i].table);
431  SG_FREE(beta_cache[i].table);
432  alpha_cache[i].table=NULL;
433  beta_cache[i].table=NULL;
434  }
435 
436  SG_FREE(alpha_cache);
437  SG_FREE(beta_cache);
438  alpha_cache=NULL;
439  beta_cache=NULL;
440 #else // USE_HMMPARALLEL_STRUCTURES
441  SG_FREE(alpha_cache.table);
442  SG_FREE(beta_cache.table);
443  alpha_cache.table=NULL;
444  beta_cache.table=NULL;
445 #endif // USE_HMMPARALLEL_STRUCTURES
446 
447  SG_FREE(states_per_observation_psi);
448  states_per_observation_psi=NULL;
449  }
450 
451 #ifdef USE_LOGSUMARRAY
452 #ifdef USE_HMMPARALLEL_STRUCTURES
453  {
454  for (int32_t i=0; i<parallel->get_num_threads(); i++)
455  SG_FREE(arrayS[i]);
456  SG_FREE(arrayS);
457  } ;
458 #else //USE_HMMPARALLEL_STRUCTURES
459  SG_FREE(arrayS);
460 #endif //USE_HMMPARALLEL_STRUCTURES
461 #endif //USE_LOGSUMARRAY
462 
463  if (!reused_caches)
464  {
465 #ifdef USE_HMMPARALLEL_STRUCTURES
466  SG_FREE(path_prob_updated);
467  SG_FREE(path_prob_dimension);
468  for (int32_t i=0; i<parallel->get_num_threads(); i++)
469  SG_FREE(path[i]);
470 #endif //USE_HMMPARALLEL_STRUCTURES
471  SG_FREE(path);
472  }
473 }
474 
475 bool CHMM::train(CFeatures* data)
476 {
477  if (data)
478  {
479  if (data->get_feature_class() != C_STRING ||
480  data->get_feature_type() != F_WORD)
481  {
482  SG_ERROR("Expected features of class string type word\n")
483  }
484  set_observations((CStringFeatures<uint16_t>*) data);
485  }
486  return baum_welch_viterbi_train(BW_NORMAL);
487 }
488 
489 bool CHMM::alloc_state_dependend_arrays()
490 {
491 
492  if (!transition_matrix_a && !observation_matrix_b &&
493  !initial_state_distribution_p && !end_state_distribution_q)
494  {
495  transition_matrix_a=SG_MALLOC(float64_t, N*N);
496  observation_matrix_b=SG_MALLOC(float64_t, N*M);
497  initial_state_distribution_p=SG_MALLOC(float64_t, N);
498  end_state_distribution_q=SG_MALLOC(float64_t, N);
499  init_model_random();
500  convert_to_log();
501  }
502 
503 #ifdef USE_HMMPARALLEL_STRUCTURES
504  for (int32_t i=0; i<parallel->get_num_threads(); i++)
505  {
506  arrayN1[i]=SG_MALLOC(float64_t, N);
507  arrayN2[i]=SG_MALLOC(float64_t, N);
508  }
509 #else //USE_HMMPARALLEL_STRUCTURES
510  arrayN1=SG_MALLOC(float64_t, N);
511  arrayN2=SG_MALLOC(float64_t, N);
512 #endif //USE_HMMPARALLEL_STRUCTURES
513 
514 #ifdef LOG_SUMARRAY
515 #ifdef USE_HMMPARALLEL_STRUCTURES
516  for (int32_t i=0; i<parallel->get_num_threads(); i++)
517  arrayS[i]=SG_MALLOC(float64_t, (int32_t)(this->N/2+1));
518 #else //USE_HMMPARALLEL_STRUCTURES
519  arrayS=SG_MALLOC(float64_t, (int32_t)(this->N/2+1));
520 #endif //USE_HMMPARALLEL_STRUCTURES
521 #endif //LOG_SUMARRAY
522  transition_matrix_A=SG_MALLOC(float64_t, this->N*this->N);
523  observation_matrix_B=SG_MALLOC(float64_t, this->N*this->M);
524 
525  if (p_observations)
526  {
527 #ifdef USE_HMMPARALLEL_STRUCTURES
528  if (alpha_cache[0].table!=NULL)
529 #else //USE_HMMPARALLEL_STRUCTURES
530  if (alpha_cache.table!=NULL)
531 #endif //USE_HMMPARALLEL_STRUCTURES
532  set_observations(p_observations);
533  else
534  set_observation_nocache(p_observations);
535  SG_UNREF(p_observations);
536  }
537 
538  this->invalidate_model();
539 
540  return ((transition_matrix_A != NULL) && (observation_matrix_B != NULL) &&
541  (transition_matrix_a != NULL) && (observation_matrix_b != NULL) &&
542  (initial_state_distribution_p != NULL) &&
543  (end_state_distribution_q != NULL));
544 }
545 
546 void CHMM::free_state_dependend_arrays()
547 {
548 #ifdef USE_HMMPARALLEL_STRUCTURES
549  for (int32_t i=0; i<parallel->get_num_threads(); i++)
550  {
551  SG_FREE(arrayN1[i]);
552  SG_FREE(arrayN2[i]);
553 
554  arrayN1[i]=NULL;
555  arrayN2[i]=NULL;
556  }
557 #else
558  SG_FREE(arrayN1);
559  SG_FREE(arrayN2);
560  arrayN1=NULL;
561  arrayN2=NULL;
562 #endif
563  if (observation_matrix_b)
564  {
565  SG_FREE(transition_matrix_A);
566  SG_FREE(observation_matrix_B);
567  SG_FREE(transition_matrix_a);
568  SG_FREE(observation_matrix_b);
569  SG_FREE(initial_state_distribution_p);
570  SG_FREE(end_state_distribution_q);
571  } ;
572 
573  transition_matrix_A=NULL;
574  observation_matrix_B=NULL;
575  transition_matrix_a=NULL;
576  observation_matrix_b=NULL;
577  initial_state_distribution_p=NULL;
578  end_state_distribution_q=NULL;
579 }
580 
581 bool CHMM::initialize(Model* m, float64_t pseudo, FILE* modelfile)
582 {
583  //yes optimistic
584  bool files_ok=true;
585 
586  trans_list_forward = NULL ;
587  trans_list_forward_cnt = NULL ;
588  trans_list_forward_val = NULL ;
589  trans_list_backward = NULL ;
590  trans_list_backward_cnt = NULL ;
591  trans_list_len = 0 ;
592  mem_initialized = false ;
593 
594  this->transition_matrix_a=NULL;
595  this->observation_matrix_b=NULL;
596  this->initial_state_distribution_p=NULL;
597  this->end_state_distribution_q=NULL;
598  this->PSEUDO= pseudo;
599  this->model= m;
600  this->p_observations=NULL;
601  this->reused_caches=false;
602 
603 #ifdef USE_HMMPARALLEL_STRUCTURES
604  alpha_cache=SG_MALLOC(T_ALPHA_BETA, parallel->get_num_threads());
605  beta_cache=SG_MALLOC(T_ALPHA_BETA, parallel->get_num_threads());
606  states_per_observation_psi=SG_MALLOC(P_STATES, parallel->get_num_threads());
607 
608  for (int32_t i=0; i<parallel->get_num_threads(); i++)
609  {
610  this->alpha_cache[i].table=NULL;
611  this->beta_cache[i].table=NULL;
612  this->alpha_cache[i].dimension=0;
613  this->beta_cache[i].dimension=0;
614  this->states_per_observation_psi[i]=NULL ;
615  }
616 
617 #else // USE_HMMPARALLEL_STRUCTURES
618  this->alpha_cache.table=NULL;
619  this->beta_cache.table=NULL;
620  this->alpha_cache.dimension=0;
621  this->beta_cache.dimension=0;
622  this->states_per_observation_psi=NULL ;
623 #endif //USE_HMMPARALLEL_STRUCTURES
624 
625  if (modelfile)
626  files_ok= files_ok && load_model(modelfile);
627 
628 #ifdef USE_HMMPARALLEL_STRUCTURES
629  path_prob_updated=SG_MALLOC(bool, parallel->get_num_threads());
630  path_prob_dimension=SG_MALLOC(int, parallel->get_num_threads());
631 
632  path=SG_MALLOC(P_STATES, parallel->get_num_threads());
633 
634  for (int32_t i=0; i<parallel->get_num_threads(); i++)
635  this->path[i]=NULL;
636 
637 #else // USE_HMMPARALLEL_STRUCTURES
638  this->path=NULL;
639 
640 #endif //USE_HMMPARALLEL_STRUCTURES
641 
642 #ifdef USE_HMMPARALLEL_STRUCTURES
643  arrayN1=SG_MALLOC(float64_t*, parallel->get_num_threads());
644  arrayN2=SG_MALLOC(float64_t*, parallel->get_num_threads());
645 #endif //USE_HMMPARALLEL_STRUCTURES
646 
647 #ifdef LOG_SUMARRAY
648 #ifdef USE_HMMPARALLEL_STRUCTURES
649  arrayS=SG_MALLOC(float64_t*, parallel->get_num_threads());
650 #endif // USE_HMMPARALLEL_STRUCTURES
651 #endif //LOG_SUMARRAY
652 
653  alloc_state_dependend_arrays();
654 
655  this->loglikelihood=false;
656  mem_initialized = true ;
657  this->invalidate_model();
658 
659  return ((files_ok) &&
660  (transition_matrix_A != NULL) && (observation_matrix_B != NULL) &&
661  (transition_matrix_a != NULL) && (observation_matrix_b != NULL) && (initial_state_distribution_p != NULL) &&
662  (end_state_distribution_q != NULL));
663 }
664 
665 //------------------------------------------------------------------------------------//
666 
667 //forward algorithm
668 //calculates Pr[O_0,O_1, ..., O_t, q_time=S_i| lambda] for 0<= time <= T-1
669 //Pr[O|lambda] for time > T
670 float64_t CHMM::forward_comp(int32_t time, int32_t state, int32_t dimension)
671 {
672  T_ALPHA_BETA_TABLE* alpha_new;
673  T_ALPHA_BETA_TABLE* alpha;
674  T_ALPHA_BETA_TABLE* dummy;
675  if (time<1)
676  time=0;
677 
678 
679  int32_t wanted_time=time;
680 
681  if (ALPHA_CACHE(dimension).table)
682  {
683  alpha=&ALPHA_CACHE(dimension).table[0];
684  alpha_new=&ALPHA_CACHE(dimension).table[N];
685  time=p_observations->get_vector_length(dimension)+1;
686  }
687  else
688  {
689  alpha_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
690  alpha=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
691  }
692 
693  if (time<1)
694  return get_p(state) + get_b(state, p_observations->get_feature(dimension,0));
695  else
696  {
697  //initialization alpha_1(i)=p_i*b_i(O_1)
698  for (int32_t i=0; i<N; i++)
699  alpha[i] = get_p(i) + get_b(i, p_observations->get_feature(dimension,0)) ;
700 
701  //induction alpha_t+1(j) = (sum_i=1^N alpha_t(i)a_ij) b_j(O_t+1)
702  for (register int32_t t=1; t<time && t < p_observations->get_vector_length(dimension); t++)
703  {
704 
705  for (int32_t j=0; j<N; j++)
706  {
707  register int32_t i, num = trans_list_forward_cnt[j] ;
708  float64_t sum=-CMath::INFTY;
709  for (i=0; i < num; i++)
710  {
711  int32_t ii = trans_list_forward[j][i] ;
712  sum = CMath::logarithmic_sum(sum, alpha[ii] + get_a(ii,j));
713  } ;
714 
715  alpha_new[j]= sum + get_b(j, p_observations->get_feature(dimension,t));
716  }
717 
718  if (!ALPHA_CACHE(dimension).table)
719  {
720  dummy=alpha;
721  alpha=alpha_new;
722  alpha_new=dummy; //switch alpha/alpha_new
723  }
724  else
725  {
726  alpha=alpha_new;
727  alpha_new+=N; //perversely pointer arithmetic
728  }
729  }
730 
731 
732  if (time<p_observations->get_vector_length(dimension))
733  {
734  register int32_t i, num=trans_list_forward_cnt[state];
735  register float64_t sum=-CMath::INFTY;
736  for (i=0; i<num; i++)
737  {
738  int32_t ii = trans_list_forward[state][i] ;
739  sum= CMath::logarithmic_sum(sum, alpha[ii] + get_a(ii, state));
740  } ;
741 
742  return sum + get_b(state, p_observations->get_feature(dimension,time));
743  }
744  else
745  {
746  // termination
747  register int32_t i ;
748  float64_t sum ;
749  sum=-CMath::INFTY;
750  for (i=0; i<N; i++) //sum over all paths
751  sum=CMath::logarithmic_sum(sum, alpha[i] + get_q(i)); //to get model probability
752 
753  if (!ALPHA_CACHE(dimension).table)
754  return sum;
755  else
756  {
757  ALPHA_CACHE(dimension).dimension=dimension;
758  ALPHA_CACHE(dimension).updated=true;
759  ALPHA_CACHE(dimension).sum=sum;
760 
761  if (wanted_time<p_observations->get_vector_length(dimension))
762  return ALPHA_CACHE(dimension).table[wanted_time*N+state];
763  else
764  return ALPHA_CACHE(dimension).sum;
765  }
766  }
767  }
768 }
769 
770 
771 //forward algorithm
772 //calculates Pr[O_0,O_1, ..., O_t, q_time=S_i| lambda] for 0<= time <= T-1
773 //Pr[O|lambda] for time > T
774 float64_t CHMM::forward_comp_old(int32_t time, int32_t state, int32_t dimension)
775 {
776  T_ALPHA_BETA_TABLE* alpha_new;
777  T_ALPHA_BETA_TABLE* alpha;
778  T_ALPHA_BETA_TABLE* dummy;
779  if (time<1)
780  time=0;
781 
782  int32_t wanted_time=time;
783 
784  if (ALPHA_CACHE(dimension).table)
785  {
786  alpha=&ALPHA_CACHE(dimension).table[0];
787  alpha_new=&ALPHA_CACHE(dimension).table[N];
788  time=p_observations->get_vector_length(dimension)+1;
789  }
790  else
791  {
792  alpha_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
793  alpha=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
794  }
795 
796  if (time<1)
797  return get_p(state) + get_b(state, p_observations->get_feature(dimension,0));
798  else
799  {
800  //initialization alpha_1(i)=p_i*b_i(O_1)
801  for (int32_t i=0; i<N; i++)
802  alpha[i] = get_p(i) + get_b(i, p_observations->get_feature(dimension,0)) ;
803 
804  //induction alpha_t+1(j) = (sum_i=1^N alpha_t(i)a_ij) b_j(O_t+1)
805  for (register int32_t t=1; t<time && t < p_observations->get_vector_length(dimension); t++)
806  {
807 
808  for (int32_t j=0; j<N; j++)
809  {
810  register int32_t i ;
811 #ifdef USE_LOGSUMARRAY
812  for (i=0; i<(N>>1); i++)
813  ARRAYS(dimension)[i]=CMath::logarithmic_sum(alpha[i<<1] + get_a(i<<1,j),
814  alpha[(i<<1)+1] + get_a((i<<1)+1,j));
815  if (N%2==1)
816  alpha_new[j]=get_b(j, p_observations->get_feature(dimension,t))+
817  CMath::logarithmic_sum(alpha[N-1]+get_a(N-1,j),
818  CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
819  else
820  alpha_new[j]=get_b(j, p_observations->get_feature(dimension,t))+CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
821 #else //USE_LOGSUMARRAY
822  float64_t sum=-CMath::INFTY;
823  for (i=0; i<N; i++)
824  sum= CMath::logarithmic_sum(sum, alpha[i] + get_a(i,j));
825 
826  alpha_new[j]= sum + get_b(j, p_observations->get_feature(dimension,t));
827 #endif //USE_LOGSUMARRAY
828  }
829 
830  if (!ALPHA_CACHE(dimension).table)
831  {
832  dummy=alpha;
833  alpha=alpha_new;
834  alpha_new=dummy; //switch alpha/alpha_new
835  }
836  else
837  {
838  alpha=alpha_new;
839  alpha_new+=N; //perversely pointer arithmetic
840  }
841  }
842 
843 
844  if (time<p_observations->get_vector_length(dimension))
845  {
846  register int32_t i;
847 #ifdef USE_LOGSUMARRAY
848  for (i=0; i<(N>>1); i++)
849  ARRAYS(dimension)[i]=CMath::logarithmic_sum(alpha[i<<1] + get_a(i<<1,state),
850  alpha[(i<<1)+1] + get_a((i<<1)+1,state));
851  if (N%2==1)
852  return get_b(state, p_observations->get_feature(dimension,time))+
853  CMath::logarithmic_sum(alpha[N-1]+get_a(N-1,state),
854  CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
855  else
856  return get_b(state, p_observations->get_feature(dimension,time))+CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
857 #else //USE_LOGSUMARRAY
858  register float64_t sum=-CMath::INFTY;
859  for (i=0; i<N; i++)
860  sum= CMath::logarithmic_sum(sum, alpha[i] + get_a(i, state));
861 
862  return sum + get_b(state, p_observations->get_feature(dimension,time));
863 #endif //USE_LOGSUMARRAY
864  }
865  else
866  {
867  // termination
868  register int32_t i ;
869  float64_t sum ;
870 #ifdef USE_LOGSUMARRAY
871  for (i=0; i<(N>>1); i++)
872  ARRAYS(dimension)[i]=CMath::logarithmic_sum(alpha[i<<1] + get_q(i<<1),
873  alpha[(i<<1)+1] + get_q((i<<1)+1));
874  if (N%2==1)
875  sum=CMath::logarithmic_sum(alpha[N-1]+get_q(N-1),
876  CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
877  else
878  sum=CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
879 #else //USE_LOGSUMARRAY
880  sum=-CMath::INFTY;
881  for (i=0; i<N; i++) //sum over all paths
882  sum=CMath::logarithmic_sum(sum, alpha[i] + get_q(i)); //to get model probability
883 #endif //USE_LOGSUMARRAY
884 
885  if (!ALPHA_CACHE(dimension).table)
886  return sum;
887  else
888  {
889  ALPHA_CACHE(dimension).dimension=dimension;
890  ALPHA_CACHE(dimension).updated=true;
891  ALPHA_CACHE(dimension).sum=sum;
892 
893  if (wanted_time<p_observations->get_vector_length(dimension))
894  return ALPHA_CACHE(dimension).table[wanted_time*N+state];
895  else
896  return ALPHA_CACHE(dimension).sum;
897  }
898  }
899  }
900 }
901 
902 
903 //backward algorithm
904 //calculates Pr[O_t+1,O_t+2, ..., O_T| q_time=S_i, lambda] for 0<= time <= T-1
905 //Pr[O|lambda] for time >= T
906 float64_t CHMM::backward_comp(int32_t time, int32_t state, int32_t dimension)
907 {
908  T_ALPHA_BETA_TABLE* beta_new;
909  T_ALPHA_BETA_TABLE* beta;
910  T_ALPHA_BETA_TABLE* dummy;
911  int32_t wanted_time=time;
912 
913  if (time<0)
914  forward(time, state, dimension);
915 
916  if (BETA_CACHE(dimension).table)
917  {
918  beta=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-1)];
919  beta_new=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-2)];
920  time=-1;
921  }
922  else
923  {
924  beta_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
925  beta=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
926  }
927 
928  if (time>=p_observations->get_vector_length(dimension)-1)
929  // return 0;
930  // else if (time==p_observations->get_vector_length(dimension)-1)
931  return get_q(state);
932  else
933  {
934  //initialization beta_T(i)=q(i)
935  for (register int32_t i=0; i<N; i++)
936  beta[i]=get_q(i);
937 
938  //induction beta_t(i) = (sum_j=1^N a_ij*b_j(O_t+1)*beta_t+1(j)
939  for (register int32_t t=p_observations->get_vector_length(dimension)-1; t>time+1 && t>0; t--)
940  {
941  for (register int32_t i=0; i<N; i++)
942  {
943  register int32_t j, num=trans_list_backward_cnt[i] ;
944  float64_t sum=-CMath::INFTY;
945  for (j=0; j<num; j++)
946  {
947  int32_t jj = trans_list_backward[i][j] ;
948  sum= CMath::logarithmic_sum(sum, get_a(i, jj) + get_b(jj, p_observations->get_feature(dimension,t)) + beta[jj]);
949  } ;
950  beta_new[i]=sum;
951  }
952 
953  if (!BETA_CACHE(dimension).table)
954  {
955  dummy=beta;
956  beta=beta_new;
957  beta_new=dummy; //switch beta/beta_new
958  }
959  else
960  {
961  beta=beta_new;
962  beta_new-=N; //perversely pointer arithmetic
963  }
964  }
965 
966  if (time>=0)
967  {
968  register int32_t j, num=trans_list_backward_cnt[state] ;
969  float64_t sum=-CMath::INFTY;
970  for (j=0; j<num; j++)
971  {
972  int32_t jj = trans_list_backward[state][j] ;
973  sum= CMath::logarithmic_sum(sum, get_a(state, jj) + get_b(jj, p_observations->get_feature(dimension,time+1))+beta[jj]);
974  } ;
975  return sum;
976  }
977  else // time<0
978  {
979  if (BETA_CACHE(dimension).table)
980  {
981  float64_t sum=-CMath::INFTY;
982  for (register int32_t j=0; j<N; j++)
983  sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
984  BETA_CACHE(dimension).sum=sum;
985  BETA_CACHE(dimension).dimension=dimension;
986  BETA_CACHE(dimension).updated=true;
987 
988  if (wanted_time<p_observations->get_vector_length(dimension))
989  return BETA_CACHE(dimension).table[wanted_time*N+state];
990  else
991  return BETA_CACHE(dimension).sum;
992  }
993  else
994  {
995  float64_t sum=-CMath::INFTY; // apply LOG_SUM_ARRAY -- no cache ... does not make very sense anyway...
996  for (register int32_t j=0; j<N; j++)
997  sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
998  return sum;
999  }
1000  }
1001  }
1002 }
1003 
1004 
1005 float64_t CHMM::backward_comp_old(
1006  int32_t time, int32_t state, int32_t dimension)
1007 {
1008  T_ALPHA_BETA_TABLE* beta_new;
1009  T_ALPHA_BETA_TABLE* beta;
1010  T_ALPHA_BETA_TABLE* dummy;
1011  int32_t wanted_time=time;
1012 
1013  if (time<0)
1014  forward(time, state, dimension);
1015 
1016  if (BETA_CACHE(dimension).table)
1017  {
1018  beta=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-1)];
1019  beta_new=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-2)];
1020  time=-1;
1021  }
1022  else
1023  {
1024  beta_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
1025  beta=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
1026  }
1027 
1028  if (time>=p_observations->get_vector_length(dimension)-1)
1029  // return 0;
1030  // else if (time==p_observations->get_vector_length(dimension)-1)
1031  return get_q(state);
1032  else
1033  {
1034  //initialization beta_T(i)=q(i)
1035  for (register int32_t i=0; i<N; i++)
1036  beta[i]=get_q(i);
1037 
1038  //induction beta_t(i) = (sum_j=1^N a_ij*b_j(O_t+1)*beta_t+1(j)
1039  for (register int32_t t=p_observations->get_vector_length(dimension)-1; t>time+1 && t>0; t--)
1040  {
1041  for (register int32_t i=0; i<N; i++)
1042  {
1043  register int32_t j ;
1044 #ifdef USE_LOGSUMARRAY
1045  for (j=0; j<(N>>1); j++)
1046  ARRAYS(dimension)[j]=CMath::logarithmic_sum(
1047  get_a(i, j<<1) + get_b(j<<1, p_observations->get_feature(dimension,t)) + beta[j<<1],
1048  get_a(i, (j<<1)+1) + get_b((j<<1)+1, p_observations->get_feature(dimension,t)) + beta[(j<<1)+1]);
1049  if (N%2==1)
1050  beta_new[i]=CMath::logarithmic_sum(get_a(i, N-1) + get_b(N-1, p_observations->get_feature(dimension,t)) + beta[N-1],
1051  CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
1052  else
1053  beta_new[i]=CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
1054 #else //USE_LOGSUMARRAY
1055  float64_t sum=-CMath::INFTY;
1056  for (j=0; j<N; j++)
1057  sum= CMath::logarithmic_sum(sum, get_a(i, j) + get_b(j, p_observations->get_feature(dimension,t)) + beta[j]);
1058 
1059  beta_new[i]=sum;
1060 #endif //USE_LOGSUMARRAY
1061  }
1062 
1063  if (!BETA_CACHE(dimension).table)
1064  {
1065  dummy=beta;
1066  beta=beta_new;
1067  beta_new=dummy; //switch beta/beta_new
1068  }
1069  else
1070  {
1071  beta=beta_new;
1072  beta_new-=N; //perversely pointer arithmetic
1073  }
1074  }
1075 
1076  if (time>=0)
1077  {
1078  register int32_t j ;
1079 #ifdef USE_LOGSUMARRAY
1080  for (j=0; j<(N>>1); j++)
1081  ARRAYS(dimension)[j]=CMath::logarithmic_sum(
1082  get_a(state, j<<1) + get_b(j<<1, p_observations->get_feature(dimension,time+1)) + beta[j<<1],
1083  get_a(state, (j<<1)+1) + get_b((j<<1)+1, p_observations->get_feature(dimension,time+1)) + beta[(j<<1)+1]);
1084  if (N%2==1)
1085  return CMath::logarithmic_sum(get_a(state, N-1) + get_b(N-1, p_observations->get_feature(dimension,time+1)) + beta[N-1],
1086  CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
1087  else
1088  return CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
1089 #else //USE_LOGSUMARRAY
1090  float64_t sum=-CMath::INFTY;
1091  for (j=0; j<N; j++)
1092  sum= CMath::logarithmic_sum(sum, get_a(state, j) + get_b(j, p_observations->get_feature(dimension,time+1))+beta[j]);
1093 
1094  return sum;
1095 #endif //USE_LOGSUMARRAY
1096  }
1097  else // time<0
1098  {
1099  if (BETA_CACHE(dimension).table)
1100  {
1101 #ifdef USE_LOGSUMARRAY//AAA
1102  for (int32_t j=0; j<(N>>1); j++)
1103  ARRAYS(dimension)[j]=CMath::logarithmic_sum(get_p(j<<1) + get_b(j<<1, p_observations->get_feature(dimension,0))+beta[j<<1],
1104  get_p((j<<1)+1) + get_b((j<<1)+1, p_observations->get_feature(dimension,0))+beta[(j<<1)+1]) ;
1105  if (N%2==1)
1106  BETA_CACHE(dimension).sum=CMath::logarithmic_sum(get_p(N-1) + get_b(N-1, p_observations->get_feature(dimension,0))+beta[N-1],
1107  CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
1108  else
1109  BETA_CACHE(dimension).sum=CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
1110 #else //USE_LOGSUMARRAY
1111  float64_t sum=-CMath::INFTY;
1112  for (register int32_t j=0; j<N; j++)
1113  sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
1114  BETA_CACHE(dimension).sum=sum;
1115 #endif //USE_LOGSUMARRAY
1116  BETA_CACHE(dimension).dimension=dimension;
1117  BETA_CACHE(dimension).updated=true;
1118 
1119  if (wanted_time<p_observations->get_vector_length(dimension))
1120  return BETA_CACHE(dimension).table[wanted_time*N+state];
1121  else
1122  return BETA_CACHE(dimension).sum;
1123  }
1124  else
1125  {
1126  float64_t sum=-CMath::INFTY; // apply LOG_SUM_ARRAY -- no cache ... does not make very sense anyway...
1127  for (register int32_t j=0; j<N; j++)
1128  sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
1129  return sum;
1130  }
1131  }
1132  }
1133 }
1134 
1135 //calculates probability of best path through the model lambda AND path itself
1136 //using viterbi algorithm
1137 float64_t CHMM::best_path(int32_t dimension)
1138 {
1139  if (!p_observations)
1140  return -1;
1141 
1142  if (dimension==-1)
1143  {
1144  if (!all_path_prob_updated)
1145  {
1146  SG_INFO("computing full viterbi likelihood\n")
1147  float64_t sum = 0 ;
1148  for (int32_t i=0; i<p_observations->get_num_vectors(); i++)
1149  sum+=best_path(i) ;
1150  sum /= p_observations->get_num_vectors() ;
1151  all_pat_prob=sum ;
1152  all_path_prob_updated=true ;
1153  return sum ;
1154  } else
1155  return all_pat_prob ;
1156  } ;
1157 
1158  if (!STATES_PER_OBSERVATION_PSI(dimension))
1159  return -1 ;
1160 
1161  if (dimension >= p_observations->get_num_vectors())
1162  return -1;
1163 
1164  if (PATH_PROB_UPDATED(dimension) && dimension==PATH_PROB_DIMENSION(dimension))
1165  return pat_prob;
1166  else
1167  {
1168  register float64_t* delta= ARRAYN2(dimension);
1169  register float64_t* delta_new= ARRAYN1(dimension);
1170 
1171  { //initialization
1172  for (register int32_t i=0; i<N; i++)
1173  {
1174  delta[i]=get_p(i)+get_b(i, p_observations->get_feature(dimension,0));
1175  set_psi(0, i, 0, dimension);
1176  }
1177  }
1178 
1179 #ifdef USE_PATHDEBUG
1180  float64_t worst=-CMath::INFTY/4 ;
1181 #endif
1182  //recursion
1183  for (register int32_t t=1; t<p_observations->get_vector_length(dimension); t++)
1184  {
1185  register float64_t* dummy;
1186  register int32_t NN=N ;
1187  for (register int32_t j=0; j<NN; j++)
1188  {
1189  register float64_t * matrix_a=&transition_matrix_a[j*N] ; // sorry for that
1190  register float64_t maxj=delta[0] + matrix_a[0];
1191  register int32_t argmax=0;
1192 
1193  for (register int32_t i=1; i<NN; i++)
1194  {
1195  register float64_t temp = delta[i] + matrix_a[i];
1196 
1197  if (temp>maxj)
1198  {
1199  maxj=temp;
1200  argmax=i;
1201  }
1202  }
1203 #ifdef FIX_POS
1204  if ((!model) || (model->get_fix_pos_state(t,j,NN)!=Model::FIX_DISALLOWED))
1205 #endif
1206  delta_new[j]=maxj + get_b(j,p_observations->get_feature(dimension,t));
1207 #ifdef FIX_POS
1208  else
1209  delta_new[j]=maxj + get_b(j,p_observations->get_feature(dimension,t)) + Model::DISALLOWED_PENALTY;
1210 #endif
1211  set_psi(t, j, argmax, dimension);
1212  }
1213 
1214 #ifdef USE_PATHDEBUG
1215  float64_t best=log(0) ;
1216  for (int32_t jj=0; jj<N; jj++)
1217  if (delta_new[jj]>best)
1218  best=delta_new[jj] ;
1219 
1220  if (best<-CMath::INFTY/2)
1221  {
1222  SG_DEBUG("worst case at %i: %e:%e\n", t, best, worst)
1223  worst=best ;
1224  } ;
1225 #endif
1226 
1227  dummy=delta;
1228  delta=delta_new;
1229  delta_new=dummy; //switch delta/delta_new
1230  }
1231 
1232 
1233  { //termination
1234  register float64_t maxj=delta[0]+get_q(0);
1235  register int32_t argmax=0;
1236 
1237  for (register int32_t i=1; i<N; i++)
1238  {
1239  register float64_t temp=delta[i]+get_q(i);
1240 
1241  if (temp>maxj)
1242  {
1243  maxj=temp;
1244  argmax=i;
1245  }
1246  }
1247  pat_prob=maxj;
1248  PATH(dimension)[p_observations->get_vector_length(dimension)-1]=argmax;
1249  } ;
1250 
1251 
1252  { //state sequence backtracking
1253  for (register int32_t t=p_observations->get_vector_length(dimension)-1; t>0; t--)
1254  {
1255  PATH(dimension)[t-1]=get_psi(t, PATH(dimension)[t], dimension);
1256  }
1257  }
1258  PATH_PROB_UPDATED(dimension)=true;
1259  PATH_PROB_DIMENSION(dimension)=dimension;
1260  return pat_prob ;
1261  }
1262 }
1263 
1264 #ifndef USE_HMMPARALLEL
1265 float64_t CHMM::model_probability_comp()
1266 {
1267  //for faster calculation cache model probability
1268  mod_prob=0 ;
1269  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++) //sum in log space
1270  mod_prob+=forward(p_observations->get_vector_length(dim), 0, dim);
1271 
1272  mod_prob_updated=true;
1273  return mod_prob;
1274 }
1275 
1276 #else
1277 
1278 float64_t CHMM::model_probability_comp()
1279 {
1280  pthread_t *threads=SG_MALLOC(pthread_t, parallel->get_num_threads());
1281  S_BW_THREAD_PARAM *params=SG_MALLOC(S_BW_THREAD_PARAM, parallel->get_num_threads());
1282 
1283  SG_INFO("computing full model probablity\n")
1284  mod_prob=0;
1285 
1286  for (int32_t cpu=0; cpu<parallel->get_num_threads(); cpu++)
1287  {
1288  params[cpu].hmm=this ;
1289  params[cpu].dim_start= p_observations->get_num_vectors()*cpu/parallel->get_num_threads();
1290  params[cpu].dim_stop= p_observations->get_num_vectors()*(cpu+1)/parallel->get_num_threads();
1291  params[cpu].p_buf=SG_MALLOC(float64_t, N);
1292  params[cpu].q_buf=SG_MALLOC(float64_t, N);
1293  params[cpu].a_buf=SG_MALLOC(float64_t, N*N);
1294  params[cpu].b_buf=SG_MALLOC(float64_t, N*M);
1295  pthread_create(&threads[cpu], NULL, bw_dim_prefetch, (void*)&params[cpu]);
1296  }
1297 
1298  for (int32_t cpu=0; cpu<parallel->get_num_threads(); cpu++)
1299  {
1300  pthread_join(threads[cpu], NULL);
1301  mod_prob+=params[cpu].ret;
1302  }
1303 
1304  for (int32_t i=0; i<parallel->get_num_threads(); i++)
1305  {
1306  SG_FREE(params[i].p_buf);
1307  SG_FREE(params[i].q_buf);
1308  SG_FREE(params[i].a_buf);
1309  SG_FREE(params[i].b_buf);
1310  }
1311 
1312  SG_FREE(threads);
1313  SG_FREE(params);
1314 
1315  mod_prob_updated=true;
1316  return mod_prob;
1317 }
1318 
1319 void* CHMM::bw_dim_prefetch(void* params)
1320 {
1321  CHMM* hmm=((S_BW_THREAD_PARAM*) params)->hmm;
1322  int32_t start=((S_BW_THREAD_PARAM*) params)->dim_start;
1323  int32_t stop=((S_BW_THREAD_PARAM*) params)->dim_stop;
1324  float64_t* p_buf=((S_BW_THREAD_PARAM*) params)->p_buf;
1325  float64_t* q_buf=((S_BW_THREAD_PARAM*) params)->q_buf;
1326  float64_t* a_buf=((S_BW_THREAD_PARAM*) params)->a_buf;
1327  float64_t* b_buf=((S_BW_THREAD_PARAM*) params)->b_buf;
1328  ((S_BW_THREAD_PARAM*)params)->ret=0;
1329 
1330  for (int32_t dim=start; dim<stop; dim++)
1331  {
1332  hmm->forward_comp(hmm->p_observations->get_vector_length(dim), hmm->N-1, dim) ;
1333  hmm->backward_comp(hmm->p_observations->get_vector_length(dim), hmm->N-1, dim) ;
1334  float64_t modprob=hmm->model_probability(dim) ;
1335  hmm->ab_buf_comp(p_buf, q_buf, a_buf, b_buf, dim) ;
1336  ((S_BW_THREAD_PARAM*)params)->ret+= modprob;
1337  }
1338  return NULL ;
1339 }
1340 
1341 void* CHMM::bw_single_dim_prefetch(void * params)
1342 {
1343  CHMM* hmm=((S_BW_THREAD_PARAM*)params)->hmm ;
1344  int32_t dim=((S_DIM_THREAD_PARAM*)params)->dim ;
1345  ((S_DIM_THREAD_PARAM*)params)->prob_sum = hmm->model_probability(dim);
1346  return NULL ;
1347 }
1348 
1349 void* CHMM::vit_dim_prefetch(void * params)
1350 {
1351  CHMM* hmm=((S_DIM_THREAD_PARAM*)params)->hmm ;
1352  int32_t dim=((S_DIM_THREAD_PARAM*)params)->dim ;
1353  ((S_DIM_THREAD_PARAM*)params)->prob_sum = hmm->best_path(dim);
1354  return NULL ;
1355 }
1356 
1357 #endif //USE_HMMPARALLEL
1358 
1359 
1360 #ifdef USE_HMMPARALLEL
1361 
1362 void CHMM::ab_buf_comp(
1363  float64_t* p_buf, float64_t* q_buf, float64_t *a_buf, float64_t* b_buf,
1364  int32_t dim)
1365 {
1366  int32_t i,j,t ;
1367  float64_t a_sum;
1368  float64_t b_sum;
1369 
1370  float64_t dimmodprob=model_probability(dim);
1371 
1372  for (i=0; i<N; i++)
1373  {
1374  //estimate initial+end state distribution numerator
1375  p_buf[i]=get_p(i)+get_b(i,p_observations->get_feature(dim,0))+backward(0,i,dim) - dimmodprob;
1376  q_buf[i]=forward(p_observations->get_vector_length(dim)-1, i, dim)+get_q(i) - dimmodprob;
1377 
1378  //estimate a
1379  for (j=0; j<N; j++)
1380  {
1381  a_sum=-CMath::INFTY;
1382 
1383  for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
1384  {
1385  a_sum= CMath::logarithmic_sum(a_sum, forward(t,i,dim)+
1386  get_a(i,j)+get_b(j,p_observations->get_feature(dim,t+1))+backward(t+1,j,dim));
1387  }
1388  a_buf[N*i+j]=a_sum-dimmodprob ;
1389  }
1390 
1391  //estimate b
1392  for (j=0; j<M; j++)
1393  {
1394  b_sum=-CMath::INFTY;
1395 
1396  for (t=0; t<p_observations->get_vector_length(dim); t++)
1397  {
1398  if (p_observations->get_feature(dim,t)==j)
1399  b_sum=CMath::logarithmic_sum(b_sum, forward(t,i,dim)+backward(t, i, dim));
1400  }
1401 
1402  b_buf[M*i+j]=b_sum-dimmodprob ;
1403  }
1404  }
1405 }
1406 
1407 //estimates new model lambda out of lambda_train using baum welch algorithm
1408 void CHMM::estimate_model_baum_welch(CHMM* hmm)
1409 {
1410  int32_t i,j,cpu;
1411  float64_t fullmodprob=0; //for all dims
1412 
1413  //clear actual model a,b,p,q are used as numerator
1414  for (i=0; i<N; i++)
1415  {
1416  if (hmm->get_p(i)>CMath::ALMOST_NEG_INFTY)
1417  set_p(i,log(PSEUDO));
1418  else
1419  set_p(i,hmm->get_p(i));
1420  if (hmm->get_q(i)>CMath::ALMOST_NEG_INFTY)
1421  set_q(i,log(PSEUDO));
1422  else
1423  set_q(i,hmm->get_q(i));
1424 
1425  for (j=0; j<N; j++)
1426  if (hmm->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
1427  set_a(i,j, log(PSEUDO));
1428  else
1429  set_a(i,j,hmm->get_a(i,j));
1430  for (j=0; j<M; j++)
1431  if (hmm->get_b(i,j)>CMath::ALMOST_NEG_INFTY)
1432  set_b(i,j, log(PSEUDO));
1433  else
1434  set_b(i,j,hmm->get_b(i,j));
1435  }
1436  invalidate_model();
1437 
1438  int32_t num_threads = parallel->get_num_threads();
1439 
1440  pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
1441  S_BW_THREAD_PARAM *params=SG_MALLOC(S_BW_THREAD_PARAM, num_threads);
1442 
1443  if (p_observations->get_num_vectors()<num_threads)
1444  num_threads=p_observations->get_num_vectors();
1445 
1446  for (cpu=0; cpu<num_threads; cpu++)
1447  {
1448  params[cpu].p_buf=SG_MALLOC(float64_t, N);
1449  params[cpu].q_buf=SG_MALLOC(float64_t, N);
1450  params[cpu].a_buf=SG_MALLOC(float64_t, N*N);
1451  params[cpu].b_buf=SG_MALLOC(float64_t, N*M);
1452 
1453  params[cpu].hmm=hmm;
1454  int32_t start = p_observations->get_num_vectors()*cpu / num_threads;
1455  int32_t stop=p_observations->get_num_vectors()*(cpu+1) / num_threads;
1456 
1457  if (cpu == parallel->get_num_threads()-1)
1458  stop=p_observations->get_num_vectors();
1459 
1460  ASSERT(start<stop)
1461  params[cpu].dim_start=start;
1462  params[cpu].dim_stop=stop;
1463 
1464  pthread_create(&threads[cpu], NULL, bw_dim_prefetch, &params[cpu]);
1465  }
1466 
1467  for (cpu=0; cpu<num_threads; cpu++)
1468  {
1469  pthread_join(threads[cpu], NULL);
1470 
1471  for (i=0; i<N; i++)
1472  {
1473  //estimate initial+end state distribution numerator
1474  set_p(i, CMath::logarithmic_sum(get_p(i), params[cpu].p_buf[i]));
1475  set_q(i, CMath::logarithmic_sum(get_q(i), params[cpu].q_buf[i]));
1476 
1477  //estimate numerator for a
1478  for (j=0; j<N; j++)
1479  set_a(i,j, CMath::logarithmic_sum(get_a(i,j), params[cpu].a_buf[N*i+j]));
1480 
1481  //estimate numerator for b
1482  for (j=0; j<M; j++)
1483  set_b(i,j, CMath::logarithmic_sum(get_b(i,j), params[cpu].b_buf[M*i+j]));
1484  }
1485 
1486  fullmodprob+=params[cpu].ret;
1487 
1488  }
1489 
1490  for (cpu=0; cpu<num_threads; cpu++)
1491  {
1492  SG_FREE(params[cpu].p_buf);
1493  SG_FREE(params[cpu].q_buf);
1494  SG_FREE(params[cpu].a_buf);
1495  SG_FREE(params[cpu].b_buf);
1496  }
1497 
1498  SG_FREE(threads);
1499  SG_FREE(params);
1500 
1501  //cache hmm model probability
1502  hmm->mod_prob=fullmodprob;
1503  hmm->mod_prob_updated=true ;
1504 
1505  //new model probability is unknown
1506  normalize();
1507  invalidate_model();
1508 }
1509 
1510 #else // USE_HMMPARALLEL
1511 
1512 //estimates new model lambda out of lambda_estimate using baum welch algorithm
1513 void CHMM::estimate_model_baum_welch(CHMM* estimate)
1514 {
1515  int32_t i,j,t,dim;
1516  float64_t a_sum, b_sum; //numerator
1517  float64_t dimmodprob=0; //model probability for dim
1518  float64_t fullmodprob=0; //for all dims
1519 
1520  //clear actual model a,b,p,q are used as numerator
1521  for (i=0; i<N; i++)
1522  {
1523  if (estimate->get_p(i)>CMath::ALMOST_NEG_INFTY)
1524  set_p(i,log(PSEUDO));
1525  else
1526  set_p(i,estimate->get_p(i));
1527  if (estimate->get_q(i)>CMath::ALMOST_NEG_INFTY)
1528  set_q(i,log(PSEUDO));
1529  else
1530  set_q(i,estimate->get_q(i));
1531 
1532  for (j=0; j<N; j++)
1533  if (estimate->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
1534  set_a(i,j, log(PSEUDO));
1535  else
1536  set_a(i,j,estimate->get_a(i,j));
1537  for (j=0; j<M; j++)
1538  if (estimate->get_b(i,j)>CMath::ALMOST_NEG_INFTY)
1539  set_b(i,j, log(PSEUDO));
1540  else
1541  set_b(i,j,estimate->get_b(i,j));
1542  }
1543  invalidate_model();
1544 
1545  //change summation order to make use of alpha/beta caches
1546  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
1547  {
1548  dimmodprob=estimate->model_probability(dim);
1549  fullmodprob+=dimmodprob ;
1550 
1551  for (i=0; i<N; i++)
1552  {
1553  //estimate initial+end state distribution numerator
1554  set_p(i, CMath::logarithmic_sum(get_p(i), estimate->get_p(i)+estimate->get_b(i,p_observations->get_feature(dim,0))+estimate->backward(0,i,dim) - dimmodprob));
1555  set_q(i, CMath::logarithmic_sum(get_q(i), estimate->forward(p_observations->get_vector_length(dim)-1, i, dim)+estimate->get_q(i) - dimmodprob ));
1556 
1557  int32_t num = trans_list_backward_cnt[i] ;
1558 
1559  //estimate a
1560  for (j=0; j<num; j++)
1561  {
1562  int32_t jj = trans_list_backward[i][j] ;
1563  a_sum=-CMath::INFTY;
1564 
1565  for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
1566  {
1567  a_sum= CMath::logarithmic_sum(a_sum, estimate->forward(t,i,dim)+
1568  estimate->get_a(i,jj)+estimate->get_b(jj,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,jj,dim));
1569  }
1570  set_a(i,jj, CMath::logarithmic_sum(get_a(i,jj), a_sum-dimmodprob));
1571  }
1572 
1573  //estimate b
1574  for (j=0; j<M; j++)
1575  {
1576  b_sum=-CMath::INFTY;
1577 
1578  for (t=0; t<p_observations->get_vector_length(dim); t++)
1579  {
1580  if (p_observations->get_feature(dim,t)==j)
1581  b_sum=CMath::logarithmic_sum(b_sum, estimate->forward(t,i,dim)+estimate->backward(t, i, dim));
1582  }
1583 
1584  set_b(i,j, CMath::logarithmic_sum(get_b(i,j), b_sum-dimmodprob));
1585  }
1586  }
1587  }
1588 
1589  //cache estimate model probability
1590  estimate->mod_prob=fullmodprob;
1591  estimate->mod_prob_updated=true ;
1592 
1593  //new model probability is unknown
1594  normalize();
1595  invalidate_model();
1596 }
1597 
1598 //estimates new model lambda out of lambda_estimate using baum welch algorithm
1599 void CHMM::estimate_model_baum_welch_old(CHMM* estimate)
1600 {
1601  int32_t i,j,t,dim;
1602  float64_t a_sum, b_sum; //numerator
1603  float64_t dimmodprob=0; //model probability for dim
1604  float64_t fullmodprob=0; //for all dims
1605 
1606  //clear actual model a,b,p,q are used as numerator
1607  for (i=0; i<N; i++)
1608  {
1609  if (estimate->get_p(i)>CMath::ALMOST_NEG_INFTY)
1610  set_p(i,log(PSEUDO));
1611  else
1612  set_p(i,estimate->get_p(i));
1613  if (estimate->get_q(i)>CMath::ALMOST_NEG_INFTY)
1614  set_q(i,log(PSEUDO));
1615  else
1616  set_q(i,estimate->get_q(i));
1617 
1618  for (j=0; j<N; j++)
1619  if (estimate->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
1620  set_a(i,j, log(PSEUDO));
1621  else
1622  set_a(i,j,estimate->get_a(i,j));
1623  for (j=0; j<M; j++)
1624  if (estimate->get_b(i,j)>CMath::ALMOST_NEG_INFTY)
1625  set_b(i,j, log(PSEUDO));
1626  else
1627  set_b(i,j,estimate->get_b(i,j));
1628  }
1629  invalidate_model();
1630 
1631  //change summation order to make use of alpha/beta caches
1632  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
1633  {
1634  dimmodprob=estimate->model_probability(dim);
1635  fullmodprob+=dimmodprob ;
1636 
1637  for (i=0; i<N; i++)
1638  {
1639  //estimate initial+end state distribution numerator
1640  set_p(i, CMath::logarithmic_sum(get_p(i), estimate->get_p(i)+estimate->get_b(i,p_observations->get_feature(dim,0))+estimate->backward(0,i,dim) - dimmodprob));
1641  set_q(i, CMath::logarithmic_sum(get_q(i), estimate->forward(p_observations->get_vector_length(dim)-1, i, dim)+estimate->get_q(i) - dimmodprob ));
1642 
1643  //estimate a
1644  for (j=0; j<N; j++)
1645  {
1646  a_sum=-CMath::INFTY;
1647 
1648  for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
1649  {
1650  a_sum= CMath::logarithmic_sum(a_sum, estimate->forward(t,i,dim)+
1651  estimate->get_a(i,j)+estimate->get_b(j,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,j,dim));
1652  }
1653  set_a(i,j, CMath::logarithmic_sum(get_a(i,j), a_sum-dimmodprob));
1654  }
1655 
1656  //estimate b
1657  for (j=0; j<M; j++)
1658  {
1659  b_sum=-CMath::INFTY;
1660 
1661  for (t=0; t<p_observations->get_vector_length(dim); t++)
1662  {
1663  if (p_observations->get_feature(dim,t)==j)
1664  b_sum=CMath::logarithmic_sum(b_sum, estimate->forward(t,i,dim)+estimate->backward(t, i, dim));
1665  }
1666 
1667  set_b(i,j, CMath::logarithmic_sum(get_b(i,j), b_sum-dimmodprob));
1668  }
1669  }
1670  }
1671 
1672  //cache estimate model probability
1673  estimate->mod_prob=fullmodprob;
1674  estimate->mod_prob_updated=true ;
1675 
1676  //new model probability is unknown
1677  normalize();
1678  invalidate_model();
1679 }
1680 #endif // USE_HMMPARALLEL
1681 
1682 //estimates new model lambda out of lambda_estimate using baum welch algorithm
1683 // optimize only p, q, a but not b
1684 void CHMM::estimate_model_baum_welch_trans(CHMM* estimate)
1685 {
1686  int32_t i,j,t,dim;
1687  float64_t a_sum; //numerator
1688  float64_t dimmodprob=0; //model probability for dim
1689  float64_t fullmodprob=0; //for all dims
1690 
1691  //clear actual model a,b,p,q are used as numerator
1692  for (i=0; i<N; i++)
1693  {
1694  if (estimate->get_p(i)>CMath::ALMOST_NEG_INFTY)
1695  set_p(i,log(PSEUDO));
1696  else
1697  set_p(i,estimate->get_p(i));
1698  if (estimate->get_q(i)>CMath::ALMOST_NEG_INFTY)
1699  set_q(i,log(PSEUDO));
1700  else
1701  set_q(i,estimate->get_q(i));
1702 
1703  for (j=0; j<N; j++)
1704  if (estimate->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
1705  set_a(i,j, log(PSEUDO));
1706  else
1707  set_a(i,j,estimate->get_a(i,j));
1708  for (j=0; j<M; j++)
1709  set_b(i,j,estimate->get_b(i,j));
1710  }
1711  invalidate_model();
1712 
1713  //change summation order to make use of alpha/beta caches
1714  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
1715  {
1716  dimmodprob=estimate->model_probability(dim);
1717  fullmodprob+=dimmodprob ;
1718 
1719  for (i=0; i<N; i++)
1720  {
1721  //estimate initial+end state distribution numerator
1722  set_p(i, CMath::logarithmic_sum(get_p(i), estimate->get_p(i)+estimate->get_b(i,p_observations->get_feature(dim,0))+estimate->backward(0,i,dim) - dimmodprob));
1723  set_q(i, CMath::logarithmic_sum(get_q(i), estimate->forward(p_observations->get_vector_length(dim)-1, i, dim)+estimate->get_q(i) - dimmodprob ));
1724 
1725  int32_t num = trans_list_backward_cnt[i] ;
1726  //estimate a
1727  for (j=0; j<num; j++)
1728  {
1729  int32_t jj = trans_list_backward[i][j] ;
1730  a_sum=-CMath::INFTY;
1731 
1732  for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
1733  {
1734  a_sum= CMath::logarithmic_sum(a_sum, estimate->forward(t,i,dim)+
1735  estimate->get_a(i,jj)+estimate->get_b(jj,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,jj,dim));
1736  }
1737  set_a(i,jj, CMath::logarithmic_sum(get_a(i,jj), a_sum-dimmodprob));
1738  }
1739  }
1740  }
1741 
1742  //cache estimate model probability
1743  estimate->mod_prob=fullmodprob;
1744  estimate->mod_prob_updated=true ;
1745 
1746  //new model probability is unknown
1747  normalize();
1748  invalidate_model();
1749 }
1750 
1751 
1752 
1753 //estimates new model lambda out of lambda_estimate using baum welch algorithm
1754 void CHMM::estimate_model_baum_welch_defined(CHMM* estimate)
1755 {
1756  int32_t i,j,old_i,k,t,dim;
1757  float64_t a_sum_num, b_sum_num; //numerator
1758  float64_t a_sum_denom, b_sum_denom; //denominator
1759  float64_t dimmodprob=-CMath::INFTY; //model probability for dim
1760  float64_t fullmodprob=0; //for all dims
1761  float64_t* A=ARRAYN1(0);
1762  float64_t* B=ARRAYN2(0);
1763 
1764  //clear actual model a,b,p,q are used as numerator
1765  //A,B as denominator for a,b
1766  for (k=0; (i=model->get_learn_p(k))!=-1; k++)
1767  set_p(i,log(PSEUDO));
1768 
1769  for (k=0; (i=model->get_learn_q(k))!=-1; k++)
1770  set_q(i,log(PSEUDO));
1771 
1772  for (k=0; (i=model->get_learn_a(k,0))!=-1; k++)
1773  {
1774  j=model->get_learn_a(k,1);
1775  set_a(i,j, log(PSEUDO));
1776  }
1777 
1778  for (k=0; (i=model->get_learn_b(k,0))!=-1; k++)
1779  {
1780  j=model->get_learn_b(k,1);
1781  set_b(i,j, log(PSEUDO));
1782  }
1783 
1784  for (i=0; i<N; i++)
1785  {
1786  A[i]=log(PSEUDO);
1787  B[i]=log(PSEUDO);
1788  }
1789 
1790 #ifdef USE_HMMPARALLEL
1791  int32_t num_threads = parallel->get_num_threads();
1792  pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
1793  S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);
1794 
1795  if (p_observations->get_num_vectors()<num_threads)
1796  num_threads=p_observations->get_num_vectors();
1797 #endif
1798 
1799  //change summation order to make use of alpha/beta caches
1800  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
1801  {
1802 #ifdef USE_HMMPARALLEL
1803  if (dim%num_threads==0)
1804  {
1805  for (i=0; i<num_threads; i++)
1806  {
1807  if (dim+i<p_observations->get_num_vectors())
1808  {
1809  params[i].hmm=estimate ;
1810  params[i].dim=dim+i ;
1811  pthread_create(&threads[i], NULL, bw_single_dim_prefetch, (void*)&params[i]) ;
1812  }
1813  }
1814  for (i=0; i<num_threads; i++)
1815  {
1816  if (dim+i<p_observations->get_num_vectors())
1817  {
1818  pthread_join(threads[i], NULL);
1819  dimmodprob = params[i].prob_sum;
1820  }
1821  }
1822  }
1823 #else
1824  dimmodprob=estimate->model_probability(dim);
1825 #endif // USE_HMMPARALLEL
1826 
1827  //and denominator
1828  fullmodprob+= dimmodprob;
1829 
1830  //estimate initial+end state distribution numerator
1831  for (k=0; (i=model->get_learn_p(k))!=-1; k++)
1832  set_p(i, CMath::logarithmic_sum(get_p(i), estimate->forward(0,i,dim)+estimate->backward(0,i,dim) - dimmodprob ) );
1833 
1834  for (k=0; (i=model->get_learn_q(k))!=-1; k++)
1835  set_q(i, CMath::logarithmic_sum(get_q(i), estimate->forward(p_observations->get_vector_length(dim)-1, i, dim)+
1836  estimate->backward(p_observations->get_vector_length(dim)-1, i, dim) - dimmodprob ) );
1837 
1838  //estimate a
1839  old_i=-1;
1840  for (k=0; (i=model->get_learn_a(k,0))!=-1; k++)
1841  {
1842  //denominator is constant for j
1843  //therefore calculate it first
1844  if (old_i!=i)
1845  {
1846  old_i=i;
1847  a_sum_denom=-CMath::INFTY;
1848 
1849  for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
1850  a_sum_denom= CMath::logarithmic_sum(a_sum_denom, estimate->forward(t,i,dim)+estimate->backward(t,i,dim));
1851 
1852  A[i]= CMath::logarithmic_sum(A[i], a_sum_denom-dimmodprob);
1853  }
1854 
1855  j=model->get_learn_a(k,1);
1856  a_sum_num=-CMath::INFTY;
1857  for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
1858  {
1859  a_sum_num= CMath::logarithmic_sum(a_sum_num, estimate->forward(t,i,dim)+
1860  estimate->get_a(i,j)+estimate->get_b(j,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,j,dim));
1861  }
1862 
1863  set_a(i,j, CMath::logarithmic_sum(get_a(i,j), a_sum_num-dimmodprob));
1864  }
1865 
1866  //estimate b
1867  old_i=-1;
1868  for (k=0; (i=model->get_learn_b(k,0))!=-1; k++)
1869  {
1870 
1871  //denominator is constant for j
1872  //therefore calculate it first
1873  if (old_i!=i)
1874  {
1875  old_i=i;
1876  b_sum_denom=-CMath::INFTY;
1877 
1878  for (t=0; t<p_observations->get_vector_length(dim); t++)
1879  b_sum_denom= CMath::logarithmic_sum(b_sum_denom, estimate->forward(t,i,dim)+estimate->backward(t,i,dim));
1880 
1881  B[i]= CMath::logarithmic_sum(B[i], b_sum_denom-dimmodprob);
1882  }
1883 
1884  j=model->get_learn_b(k,1);
1885  b_sum_num=-CMath::INFTY;
1886  for (t=0; t<p_observations->get_vector_length(dim); t++)
1887  {
1888  if (p_observations->get_feature(dim,t)==j)
1889  b_sum_num=CMath::logarithmic_sum(b_sum_num, estimate->forward(t,i,dim)+estimate->backward(t, i, dim));
1890  }
1891 
1892  set_b(i,j, CMath::logarithmic_sum(get_b(i,j), b_sum_num-dimmodprob));
1893  }
1894  }
1895 #ifdef USE_HMMPARALLEL
1896  SG_FREE(threads);
1897  SG_FREE(params);
1898 #endif
1899 
1900 
1901  //calculate estimates
1902  for (k=0; (i=model->get_learn_p(k))!=-1; k++)
1903  set_p(i, get_p(i)-log(p_observations->get_num_vectors()+N*PSEUDO) );
1904 
1905  for (k=0; (i=model->get_learn_q(k))!=-1; k++)
1906  set_q(i, get_q(i)-log(p_observations->get_num_vectors()+N*PSEUDO) );
1907 
1908  for (k=0; (i=model->get_learn_a(k,0))!=-1; k++)
1909  {
1910  j=model->get_learn_a(k,1);
1911  set_a(i,j, get_a(i,j) - A[i]);
1912  }
1913 
1914  for (k=0; (i=model->get_learn_b(k,0))!=-1; k++)
1915  {
1916  j=model->get_learn_b(k,1);
1917  set_b(i,j, get_b(i,j) - B[i]);
1918  }
1919 
1920  //cache estimate model probability
1921  estimate->mod_prob=fullmodprob;
1922  estimate->mod_prob_updated=true ;
1923 
1924  //new model probability is unknown
1925  normalize();
1926  invalidate_model();
1927 }
1928 
1929 //estimates new model lambda out of lambda_estimate using viterbi algorithm
1930 void CHMM::estimate_model_viterbi(CHMM* estimate)
1931 {
1932  int32_t i,j,t;
1933  float64_t sum;
1934  float64_t* P=ARRAYN1(0);
1935  float64_t* Q=ARRAYN2(0);
1936 
1937  path_deriv_updated=false ;
1938 
1939  //initialize with pseudocounts
1940  for (i=0; i<N; i++)
1941  {
1942  for (j=0; j<N; j++)
1943  set_A(i,j, PSEUDO);
1944 
1945  for (j=0; j<M; j++)
1946  set_B(i,j, PSEUDO);
1947 
1948  P[i]=PSEUDO;
1949  Q[i]=PSEUDO;
1950  }
1951 
1952  float64_t allpatprob=0 ;
1953 
1954 #ifdef USE_HMMPARALLEL
1955  int32_t num_threads = parallel->get_num_threads();
1956  pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
1957  S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);
1958 
1959  if (p_observations->get_num_vectors()<num_threads)
1960  num_threads=p_observations->get_num_vectors();
1961 #endif
1962 
1963  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
1964  {
1965 
1966 #ifdef USE_HMMPARALLEL
1967  if (dim%num_threads==0)
1968  {
1969  for (i=0; i<num_threads; i++)
1970  {
1971  if (dim+i<p_observations->get_num_vectors())
1972  {
1973  params[i].hmm=estimate ;
1974  params[i].dim=dim+i ;
1975  pthread_create(&threads[i], NULL, vit_dim_prefetch, (void*)&params[i]) ;
1976  }
1977  }
1978  for (i=0; i<num_threads; i++)
1979  {
1980  if (dim+i<p_observations->get_num_vectors())
1981  {
1982  pthread_join(threads[i], NULL);
1983  allpatprob += params[i].prob_sum;
1984  }
1985  }
1986  }
1987 #else
1988  //using viterbi to find best path
1989  allpatprob += estimate->best_path(dim);
1990 #endif // USE_HMMPARALLEL
1991 
1992  //counting occurences for A and B
1993  for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
1994  {
1995  set_A(estimate->PATH(dim)[t], estimate->PATH(dim)[t+1], get_A(estimate->PATH(dim)[t], estimate->PATH(dim)[t+1])+1);
1996  set_B(estimate->PATH(dim)[t], p_observations->get_feature(dim,t), get_B(estimate->PATH(dim)[t], p_observations->get_feature(dim,t))+1);
1997  }
1998 
1999  set_B(estimate->PATH(dim)[p_observations->get_vector_length(dim)-1], p_observations->get_feature(dim,p_observations->get_vector_length(dim)-1), get_B(estimate->PATH(dim)[p_observations->get_vector_length(dim)-1], p_observations->get_feature(dim,p_observations->get_vector_length(dim)-1)) + 1 );
2000 
2001  P[estimate->PATH(dim)[0]]++;
2002  Q[estimate->PATH(dim)[p_observations->get_vector_length(dim)-1]]++;
2003  }
2004 
2005 #ifdef USE_HMMPARALLEL
2006  SG_FREE(threads);
2007  SG_FREE(params);
2008 #endif
2009 
2010  allpatprob/=p_observations->get_num_vectors() ;
2011  estimate->all_pat_prob=allpatprob ;
2012  estimate->all_path_prob_updated=true ;
2013 
2014  //converting A to probability measure a
2015  for (i=0; i<N; i++)
2016  {
2017  sum=0;
2018  for (j=0; j<N; j++)
2019  sum+=get_A(i,j);
2020 
2021  for (j=0; j<N; j++)
2022  set_a(i,j, log(get_A(i,j)/sum));
2023  }
2024 
2025  //converting B to probability measures b
2026  for (i=0; i<N; i++)
2027  {
2028  sum=0;
2029  for (j=0; j<M; j++)
2030  sum+=get_B(i,j);
2031 
2032  for (j=0; j<M; j++)
2033  set_b(i,j, log(get_B(i, j)/sum));
2034  }
2035 
2036  //converting P to probability measure p
2037  sum=0;
2038  for (i=0; i<N; i++)
2039  sum+=P[i];
2040 
2041  for (i=0; i<N; i++)
2042  set_p(i, log(P[i]/sum));
2043 
2044  //converting Q to probability measure q
2045  sum=0;
2046  for (i=0; i<N; i++)
2047  sum+=Q[i];
2048 
2049  for (i=0; i<N; i++)
2050  set_q(i, log(Q[i]/sum));
2051 
2052  //new model probability is unknown
2053  invalidate_model();
2054 }
2055 
2056 // estimate parameters listed in learn_x
2057 void CHMM::estimate_model_viterbi_defined(CHMM* estimate)
2058 {
2059  int32_t i,j,k,t;
2060  float64_t sum;
2061  float64_t* P=ARRAYN1(0);
2062  float64_t* Q=ARRAYN2(0);
2063 
2064  path_deriv_updated=false ;
2065 
2066  //initialize with pseudocounts
2067  for (i=0; i<N; i++)
2068  {
2069  for (j=0; j<N; j++)
2070  set_A(i,j, PSEUDO);
2071 
2072  for (j=0; j<M; j++)
2073  set_B(i,j, PSEUDO);
2074 
2075  P[i]=PSEUDO;
2076  Q[i]=PSEUDO;
2077  }
2078 
2079 #ifdef USE_HMMPARALLEL
2080  int32_t num_threads = parallel->get_num_threads();
2081  pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
2082  S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);
2083 #endif
2084 
2085  float64_t allpatprob=0.0 ;
2086  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
2087  {
2088 
2089 #ifdef USE_HMMPARALLEL
2090  if (dim%num_threads==0)
2091  {
2092  for (i=0; i<num_threads; i++)
2093  {
2094  if (dim+i<p_observations->get_num_vectors())
2095  {
2096  params[i].hmm=estimate ;
2097  params[i].dim=dim+i ;
2098  pthread_create(&threads[i], NULL, vit_dim_prefetch, (void*)&params[i]) ;
2099  }
2100  }
2101  for (i=0; i<num_threads; i++)
2102  {
2103  if (dim+i<p_observations->get_num_vectors())
2104  {
2105  pthread_join(threads[i], NULL);
2106  allpatprob += params[i].prob_sum;
2107  }
2108  }
2109  }
2110 #else // USE_HMMPARALLEL
2111  //using viterbi to find best path
2112  allpatprob += estimate->best_path(dim);
2113 #endif // USE_HMMPARALLEL
2114 
2115 
2116  //counting occurences for A and B
2117  for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
2118  {
2119  set_A(estimate->PATH(dim)[t], estimate->PATH(dim)[t+1], get_A(estimate->PATH(dim)[t], estimate->PATH(dim)[t+1])+1);
2120  set_B(estimate->PATH(dim)[t], p_observations->get_feature(dim,t), get_B(estimate->PATH(dim)[t], p_observations->get_feature(dim,t))+1);
2121  }
2122 
2123  set_B(estimate->PATH(dim)[p_observations->get_vector_length(dim)-1], p_observations->get_feature(dim,p_observations->get_vector_length(dim)-1), get_B(estimate->PATH(dim)[p_observations->get_vector_length(dim)-1], p_observations->get_feature(dim,p_observations->get_vector_length(dim)-1)) + 1 );
2124 
2125  P[estimate->PATH(dim)[0]]++;
2126  Q[estimate->PATH(dim)[p_observations->get_vector_length(dim)-1]]++;
2127  }
2128 
2129 #ifdef USE_HMMPARALLEL
2130  SG_FREE(threads);
2131  SG_FREE(params);
2132 #endif
2133 
2134  //estimate->invalidate_model() ;
2135  //float64_t q=estimate->best_path(-1) ;
2136 
2137  allpatprob/=p_observations->get_num_vectors() ;
2138  estimate->all_pat_prob=allpatprob ;
2139  estimate->all_path_prob_updated=true ;
2140 
2141 
2142  //copy old model
2143  for (i=0; i<N; i++)
2144  {
2145  for (j=0; j<N; j++)
2146  set_a(i,j, estimate->get_a(i,j));
2147 
2148  for (j=0; j<M; j++)
2149  set_b(i,j, estimate->get_b(i,j));
2150  }
2151 
2152  //converting A to probability measure a
2153  i=0;
2154  sum=0;
2155  j=model->get_learn_a(i,0);
2156  k=i;
2157  while (model->get_learn_a(i,0)!=-1 || k<i)
2158  {
2159  if (j==model->get_learn_a(i,0))
2160  {
2161  sum+=get_A(model->get_learn_a(i,0), model->get_learn_a(i,1));
2162  i++;
2163  }
2164  else
2165  {
2166  while (k<i)
2167  {
2168  set_a(model->get_learn_a(k,0), model->get_learn_a(k,1), log (get_A(model->get_learn_a(k,0), model->get_learn_a(k,1)) / sum));
2169  k++;
2170  }
2171 
2172  sum=0;
2173  j=model->get_learn_a(i,0);
2174  k=i;
2175  }
2176  }
2177 
2178  //converting B to probability measures b
2179  i=0;
2180  sum=0;
2181  j=model->get_learn_b(i,0);
2182  k=i;
2183  while (model->get_learn_b(i,0)!=-1 || k<i)
2184  {
2185  if (j==model->get_learn_b(i,0))
2186  {
2187  sum+=get_B(model->get_learn_b(i,0),model->get_learn_b(i,1));
2188  i++;
2189  }
2190  else
2191  {
2192  while (k<i)
2193  {
2194  set_b(model->get_learn_b(k,0),model->get_learn_b(k,1), log (get_B(model->get_learn_b(k,0), model->get_learn_b(k,1)) / sum));
2195  k++;
2196  }
2197 
2198  sum=0;
2199  j=model->get_learn_b(i,0);
2200  k=i;
2201  }
2202  }
2203 
2204  i=0;
2205  sum=0;
2206  while (model->get_learn_p(i)!=-1)
2207  {
2208  sum+=P[model->get_learn_p(i)] ;
2209  i++ ;
2210  } ;
2211  i=0 ;
2212  while (model->get_learn_p(i)!=-1)
2213  {
2214  set_p(model->get_learn_p(i), log(P[model->get_learn_p(i)]/sum));
2215  i++ ;
2216  } ;
2217 
2218  i=0;
2219  sum=0;
2220  while (model->get_learn_q(i)!=-1)
2221  {
2222  sum+=Q[model->get_learn_q(i)] ;
2223  i++ ;
2224  } ;
2225  i=0 ;
2226  while (model->get_learn_q(i)!=-1)
2227  {
2228  set_q(model->get_learn_q(i), log(Q[model->get_learn_q(i)]/sum));
2229  i++ ;
2230  } ;
2231 
2232 
2233  //new model probability is unknown
2234  invalidate_model();
2235 }
2236 //------------------------------------------------------------------------------------//
2237 
2238 //to give an idea what the model looks like
2239 void CHMM::output_model(bool verbose)
2240 {
2241  int32_t i,j;
2242  float64_t checksum;
2243 
2244  //generic info
2245  SG_INFO("log(Pr[O|model])=%e, #states: %i, #observationssymbols: %i, #observations: %ix%i\n",
2246  (float64_t)((p_observations) ? model_probability() : -CMath::INFTY),
2247  N, M, ((p_observations) ? p_observations->get_max_vector_length() : 0), ((p_observations) ? p_observations->get_num_vectors() : 0));
2248 
2249  if (verbose)
2250  {
2251  // tranisition matrix a
2252  SG_INFO("\ntransition matrix\n")
2253  for (i=0; i<N; i++)
2254  {
2255  checksum= get_q(i);
2256  for (j=0; j<N; j++)
2257  {
2258  checksum= CMath::logarithmic_sum(checksum, get_a(i,j));
2259 
2260  SG_INFO("a(%02i,%02i)=%1.4f ",i,j, (float32_t) exp(get_a(i,j)))
2261 
2262  if (j % 4 == 3)
2263  SG_PRINT("\n")
2264  }
2265  if (fabs(checksum)>1e-5)
2266  SG_DEBUG(" checksum % E ******* \n",checksum)
2267  else
2268  SG_DEBUG(" checksum % E\n",checksum)
2269  }
2270 
2271  // distribution of start states p
2272  SG_INFO("\ndistribution of start states\n")
2273  checksum=-CMath::INFTY;
2274  for (i=0; i<N; i++)
2275  {
2276  checksum= CMath::logarithmic_sum(checksum, get_p(i));
2277  SG_INFO("p(%02i)=%1.4f ",i, (float32_t) exp(get_p(i)))
2278  if (i % 4 == 3)
2279  SG_PRINT("\n")
2280  }
2281  if (fabs(checksum)>1e-5)
2282  SG_DEBUG(" checksum % E ******* \n",checksum)
2283  else
2284  SG_DEBUG(" checksum=% E\n", checksum)
2285 
2286  // distribution of terminal states p
2287  SG_INFO("\ndistribution of terminal states\n")
2288  checksum=-CMath::INFTY;
2289  for (i=0; i<N; i++)
2290  {
2291  checksum= CMath::logarithmic_sum(checksum, get_q(i));
2292  SG_INFO("q(%02i)=%1.4f ",i, (float32_t) exp(get_q(i)))
2293  if (i % 4 == 3)
2294  SG_INFO("\n")
2295  }
2296  if (fabs(checksum)>1e-5)
2297  SG_DEBUG(" checksum % E ******* \n",checksum)
2298  else
2299  SG_DEBUG(" checksum=% E\n", checksum)
2300 
2301  // distribution of observations given the state b
2302  SG_INFO("\ndistribution of observations given the state\n")
2303  for (i=0; i<N; i++)
2304  {
2305  checksum=-CMath::INFTY;
2306  for (j=0; j<M; j++)
2307  {
2308  checksum=CMath::logarithmic_sum(checksum, get_b(i,j));
2309  SG_INFO("b(%02i,%02i)=%1.4f ",i,j, (float32_t) exp(get_b(i,j)))
2310  if (j % 4 == 3)
2311  SG_PRINT("\n")
2312  }
2313  if (fabs(checksum)>1e-5)
2314  SG_DEBUG(" checksum % E ******* \n",checksum)
2315  else
2316  SG_DEBUG(" checksum % E\n",checksum)
2317  }
2318  }
2319  SG_PRINT("\n")
2320 }
2321 
2322 //to give an idea what the model looks like
2323 void CHMM::output_model_defined(bool verbose)
2324 {
2325  int32_t i,j;
2326  if (!model)
2327  return ;
2328 
2329  //generic info
2330  SG_INFO("log(Pr[O|model])=%e, #states: %i, #observationssymbols: %i, #observations: %ix%i\n",
2331  (float64_t)((p_observations) ? model_probability() : -CMath::INFTY),
2332  N, M, ((p_observations) ? p_observations->get_max_vector_length() : 0), ((p_observations) ? p_observations->get_num_vectors() : 0));
2333 
2334  if (verbose)
2335  {
2336  // tranisition matrix a
2337  SG_INFO("\ntransition matrix\n")
2338 
2339  //initialize a values that have to be learned
2340  i=0;
2341  j=model->get_learn_a(i,0);
2342  while (model->get_learn_a(i,0)!=-1)
2343  {
2344  if (j!=model->get_learn_a(i,0))
2345  {
2346  j=model->get_learn_a(i,0);
2347  SG_PRINT("\n")
2348  }
2349 
2350  SG_INFO("a(%02i,%02i)=%1.4f ",model->get_learn_a(i,0), model->get_learn_a(i,1), (float32_t) exp(get_a(model->get_learn_a(i,0), model->get_learn_a(i,1))))
2351  i++;
2352  }
2353 
2354  // distribution of observations given the state b
2355  SG_INFO("\n\ndistribution of observations given the state\n")
2356  i=0;
2357  j=model->get_learn_b(i,0);
2358  while (model->get_learn_b(i,0)!=-1)
2359  {
2360  if (j!=model->get_learn_b(i,0))
2361  {
2362  j=model->get_learn_b(i,0);
2363  SG_PRINT("\n")
2364  }
2365 
2366  SG_INFO("b(%02i,%02i)=%1.4f ",model->get_learn_b(i,0),model->get_learn_b(i,1), (float32_t) exp(get_b(model->get_learn_b(i,0),model->get_learn_b(i,1))))
2367  i++;
2368  }
2369 
2370  SG_PRINT("\n")
2371  }
2372  SG_PRINT("\n")
2373 }
2374 
2375 //------------------------------------------------------------------------------------//
2376 
2377 //convert model to log probabilities
2378 void CHMM::convert_to_log()
2379 {
2380  int32_t i,j;
2381 
2382  for (i=0; i<N; i++)
2383  {
2384  if (get_p(i)!=0)
2385  set_p(i, log(get_p(i)));
2386  else
2387  set_p(i, -CMath::INFTY);;
2388  }
2389 
2390  for (i=0; i<N; i++)
2391  {
2392  if (get_q(i)!=0)
2393  set_q(i, log(get_q(i)));
2394  else
2395  set_q(i, -CMath::INFTY);;
2396  }
2397 
2398  for (i=0; i<N; i++)
2399  {
2400  for (j=0; j<N; j++)
2401  {
2402  if (get_a(i,j)!=0)
2403  set_a(i,j, log(get_a(i,j)));
2404  else
2405  set_a(i,j, -CMath::INFTY);
2406  }
2407  }
2408 
2409  for (i=0; i<N; i++)
2410  {
2411  for (j=0; j<M; j++)
2412  {
2413  if (get_b(i,j)!=0)
2414  set_b(i,j, log(get_b(i,j)));
2415  else
2416  set_b(i,j, -CMath::INFTY);
2417  }
2418  }
2419  loglikelihood=true;
2420 
2421  invalidate_model();
2422 }
2423 
2424 //init model with random values
2425 void CHMM::init_model_random()
2426 {
2427  const float64_t MIN_RAND=23e-3;
2428 
2429  float64_t sum;
2430  int32_t i,j;
2431 
2432  //initialize a with random values
2433  for (i=0; i<N; i++)
2434  {
2435  sum=0;
2436  for (j=0; j<N; j++)
2437  {
2438  set_a(i,j, CMath::random(MIN_RAND, 1.0));
2439 
2440  sum+=get_a(i,j);
2441  }
2442 
2443  for (j=0; j<N; j++)
2444  set_a(i,j, get_a(i,j)/sum);
2445  }
2446 
2447  //initialize pi with random values
2448  sum=0;
2449  for (i=0; i<N; i++)
2450  {
2451  set_p(i, CMath::random(MIN_RAND, 1.0));
2452 
2453  sum+=get_p(i);
2454  }
2455 
2456  for (i=0; i<N; i++)
2457  set_p(i, get_p(i)/sum);
2458 
2459  //initialize q with random values
2460  sum=0;
2461  for (i=0; i<N; i++)
2462  {
2463  set_q(i, CMath::random(MIN_RAND, 1.0));
2464 
2465  sum+=get_q(i);
2466  }
2467 
2468  for (i=0; i<N; i++)
2469  set_q(i, get_q(i)/sum);
2470 
2471  //initialize b with random values
2472  for (i=0; i<N; i++)
2473  {
2474  sum=0;
2475  for (j=0; j<M; j++)
2476  {
2477  set_b(i,j, CMath::random(MIN_RAND, 1.0));
2478 
2479  sum+=get_b(i,j);
2480  }
2481 
2482  for (j=0; j<M; j++)
2483  set_b(i,j, get_b(i,j)/sum);
2484  }
2485 
2486  //initialize pat/mod_prob as not calculated
2487  invalidate_model();
2488 }
2489 
2490 //init model according to const_x
2491 void CHMM::init_model_defined()
2492 {
2493  int32_t i,j,k,r;
2494  float64_t sum;
2495  const float64_t MIN_RAND=23e-3;
2496 
2497  //initialize a with zeros
2498  for (i=0; i<N; i++)
2499  for (j=0; j<N; j++)
2500  set_a(i,j, 0);
2501 
2502  //initialize p with zeros
2503  for (i=0; i<N; i++)
2504  set_p(i, 0);
2505 
2506  //initialize q with zeros
2507  for (i=0; i<N; i++)
2508  set_q(i, 0);
2509 
2510  //initialize b with zeros
2511  for (i=0; i<N; i++)
2512  for (j=0; j<M; j++)
2513  set_b(i,j, 0);
2514 
2515 
2516  //initialize a values that have to be learned
2517  float64_t *R=SG_MALLOC(float64_t, N);
2518  for (r=0; r<N; r++) R[r]=CMath::random(MIN_RAND,1.0);
2519  i=0; sum=0; k=i;
2520  j=model->get_learn_a(i,0);
2521  while (model->get_learn_a(i,0)!=-1 || k<i)
2522  {
2523  if (j==model->get_learn_a(i,0))
2524  {
2525  sum+=R[model->get_learn_a(i,1)] ;
2526  i++;
2527  }
2528  else
2529  {
2530  while (k<i)
2531  {
2532  set_a(model->get_learn_a(k,0), model->get_learn_a(k,1),
2533  R[model->get_learn_a(k,1)]/sum);
2534  k++;
2535  }
2536  j=model->get_learn_a(i,0);
2537  k=i;
2538  sum=0;
2539  for (r=0; r<N; r++) R[r]=CMath::random(MIN_RAND,1.0);
2540  }
2541  }
2542  SG_FREE(R); R=NULL ;
2543 
2544  //initialize b values that have to be learned
2545  R=SG_MALLOC(float64_t, M);
2546  for (r=0; r<M; r++) R[r]=CMath::random(MIN_RAND,1.0);
2547  i=0; sum=0; k=0 ;
2548  j=model->get_learn_b(i,0);
2549  while (model->get_learn_b(i,0)!=-1 || k<i)
2550  {
2551  if (j==model->get_learn_b(i,0))
2552  {
2553  sum+=R[model->get_learn_b(i,1)] ;
2554  i++;
2555  }
2556  else
2557  {
2558  while (k<i)
2559  {
2560  set_b(model->get_learn_b(k,0),model->get_learn_b(k,1),
2561  R[model->get_learn_b(k,1)]/sum);
2562  k++;
2563  }
2564 
2565  j=model->get_learn_b(i,0);
2566  k=i;
2567  sum=0;
2568  for (r=0; r<M; r++) R[r]=CMath::random(MIN_RAND,1.0);
2569  }
2570  }
2571  SG_FREE(R); R=NULL ;
2572 
2573  //set consts into a
2574  i=0;
2575  while (model->get_const_a(i,0) != -1)
2576  {
2577  set_a(model->get_const_a(i,0), model->get_const_a(i,1), model->get_const_a_val(i));
2578  i++;
2579  }
2580 
2581  //set consts into b
2582  i=0;
2583  while (model->get_const_b(i,0) != -1)
2584  {
2585  set_b(model->get_const_b(i,0), model->get_const_b(i,1), model->get_const_b_val(i));
2586  i++;
2587  }
2588 
2589  //set consts into p
2590  i=0;
2591  while (model->get_const_p(i) != -1)
2592  {
2593  set_p(model->get_const_p(i), model->get_const_p_val(i));
2594  i++;
2595  }
2596 
2597  //initialize q with zeros
2598  for (i=0; i<N; i++)
2599  set_q(i, 0.0);
2600 
2601  //set consts into q
2602  i=0;
2603  while (model->get_const_q(i) != -1)
2604  {
2605  set_q(model->get_const_q(i), model->get_const_q_val(i));
2606  i++;
2607  }
2608 
2609  // init p
2610  i=0;
2611  sum=0;
2612  while (model->get_learn_p(i)!=-1)
2613  {
2614  set_p(model->get_learn_p(i),CMath::random(MIN_RAND,1.0)) ;
2615  sum+=get_p(model->get_learn_p(i)) ;
2616  i++ ;
2617  } ;
2618  i=0 ;
2619  while (model->get_learn_p(i)!=-1)
2620  {
2621  set_p(model->get_learn_p(i), get_p(model->get_learn_p(i))/sum);
2622  i++ ;
2623  } ;
2624 
2625  // initialize q
2626  i=0;
2627  sum=0;
2628  while (model->get_learn_q(i)!=-1)
2629  {
2630  set_q(model->get_learn_q(i),CMath::random(MIN_RAND,1.0)) ;
2631  sum+=get_q(model->get_learn_q(i)) ;
2632  i++ ;
2633  } ;
2634  i=0 ;
2635  while (model->get_learn_q(i)!=-1)
2636  {
2637  set_q(model->get_learn_q(i), get_q(model->get_learn_q(i))/sum);
2638  i++ ;
2639  } ;
2640 
2641  //initialize pat/mod_prob as not calculated
2642  invalidate_model();
2643 }
2644 
2645 void CHMM::clear_model()
2646 {
2647  int32_t i,j;
2648  for (i=0; i<N; i++)
2649  {
2650  set_p(i, log(PSEUDO));
2651  set_q(i, log(PSEUDO));
2652 
2653  for (j=0; j<N; j++)
2654  set_a(i,j, log(PSEUDO));
2655 
2656  for (j=0; j<M; j++)
2657  set_b(i,j, log(PSEUDO));
2658  }
2659 }
2660 
2661 void CHMM::clear_model_defined()
2662 {
2663  int32_t i,j,k;
2664 
2665  for (i=0; (j=model->get_learn_p(i))!=-1; i++)
2666  set_p(j, log(PSEUDO));
2667 
2668  for (i=0; (j=model->get_learn_q(i))!=-1; i++)
2669  set_q(j, log(PSEUDO));
2670 
2671  for (i=0; (j=model->get_learn_a(i,0))!=-1; i++)
2672  {
2673  k=model->get_learn_a(i,1); // catch (j,k) as indizes to be learned
2674  set_a(j,k, log(PSEUDO));
2675  }
2676 
2677  for (i=0; (j=model->get_learn_b(i,0))!=-1; i++)
2678  {
2679  k=model->get_learn_b(i,1); // catch (j,k) as indizes to be learned
2680  set_b(j,k, log(PSEUDO));
2681  }
2682 }
2683 
2684 void CHMM::copy_model(CHMM* l)
2685 {
2686  int32_t i,j;
2687  for (i=0; i<N; i++)
2688  {
2689  set_p(i, l->get_p(i));
2690  set_q(i, l->get_q(i));
2691 
2692  for (j=0; j<N; j++)
2693  set_a(i,j, l->get_a(i,j));
2694 
2695  for (j=0; j<M; j++)
2696  set_b(i,j, l->get_b(i,j));
2697  }
2698 }
2699 
2700 void CHMM::invalidate_model()
2701 {
2702  //initialize pat/mod_prob/alpha/beta cache as not calculated
2703  this->mod_prob=0.0;
2704  this->mod_prob_updated=false;
2705 
2706  if (mem_initialized)
2707  {
2708  if (trans_list_forward_cnt)
2709  SG_FREE(trans_list_forward_cnt);
2710  trans_list_forward_cnt=NULL ;
2711  if (trans_list_backward_cnt)
2712  SG_FREE(trans_list_backward_cnt);
2713  trans_list_backward_cnt=NULL ;
2714  if (trans_list_forward)
2715  {
2716  for (int32_t i=0; i<trans_list_len; i++)
2717  if (trans_list_forward[i])
2718  SG_FREE(trans_list_forward[i]);
2719  SG_FREE(trans_list_forward);
2720  trans_list_forward=NULL ;
2721  }
2722  if (trans_list_backward)
2723  {
2724  for (int32_t i=0; i<trans_list_len; i++)
2725  if (trans_list_backward[i])
2726  SG_FREE(trans_list_backward[i]);
2727  SG_FREE(trans_list_backward);
2728  trans_list_backward = NULL ;
2729  } ;
2730 
2731  trans_list_len = N ;
2732  trans_list_forward = SG_MALLOC(T_STATES*, N);
2733  trans_list_forward_cnt = SG_MALLOC(T_STATES, N);
2734 
2735  for (int32_t j=0; j<N; j++)
2736  {
2737  trans_list_forward_cnt[j]= 0 ;
2738  trans_list_forward[j]= SG_MALLOC(T_STATES, N);
2739  for (int32_t i=0; i<N; i++)
2740  if (get_a(i,j)>CMath::ALMOST_NEG_INFTY)
2741  {
2742  trans_list_forward[j][trans_list_forward_cnt[j]]=i ;
2743  trans_list_forward_cnt[j]++ ;
2744  }
2745  } ;
2746 
2747  trans_list_backward = SG_MALLOC(T_STATES*, N);
2748  trans_list_backward_cnt = SG_MALLOC(T_STATES, N);
2749 
2750  for (int32_t i=0; i<N; i++)
2751  {
2752  trans_list_backward_cnt[i]= 0 ;
2753  trans_list_backward[i]= SG_MALLOC(T_STATES, N);
2754  for (int32_t j=0; j<N; j++)
2755  if (get_a(i,j)>CMath::ALMOST_NEG_INFTY)
2756  {
2757  trans_list_backward[i][trans_list_backward_cnt[i]]=j ;
2758  trans_list_backward_cnt[i]++ ;
2759  }
2760  } ;
2761  } ;
2762  this->all_pat_prob=0.0;
2763  this->pat_prob=0.0;
2764  this->path_deriv_updated=false ;
2765  this->path_deriv_dimension=-1 ;
2766  this->all_path_prob_updated=false;
2767 
2768 #ifdef USE_HMMPARALLEL_STRUCTURES
2769  {
2770  for (int32_t i=0; i<parallel->get_num_threads(); i++)
2771  {
2772  this->alpha_cache[i].updated=false;
2773  this->beta_cache[i].updated=false;
2774  path_prob_updated[i]=false ;
2775  path_prob_dimension[i]=-1 ;
2776  } ;
2777  }
2778 #else // USE_HMMPARALLEL_STRUCTURES
2779  this->alpha_cache.updated=false;
2780  this->beta_cache.updated=false;
2781  this->path_prob_dimension=-1;
2782  this->path_prob_updated=false;
2783 
2784 #endif // USE_HMMPARALLEL_STRUCTURES
2785 }
2786 
2787 void CHMM::open_bracket(FILE* file)
2788 {
2789  int32_t value;
2790  while (((value=fgetc(file)) != EOF) && (value!='[')) //skip possible spaces and end if '[' occurs
2791  {
2792  if (value=='\n')
2793  line++;
2794  }
2795 
2796  if (value==EOF)
2797  error(line, "expected \"[\" in input file");
2798 
2799  while (((value=fgetc(file)) != EOF) && (isspace(value))) //skip possible spaces
2800  {
2801  if (value=='\n')
2802  line++;
2803  }
2804 
2805  ungetc(value, file);
2806 }
2807 
2808 void CHMM::close_bracket(FILE* file)
2809 {
2810  int32_t value;
2811  while (((value=fgetc(file)) != EOF) && (value!=']')) //skip possible spaces and end if ']' occurs
2812  {
2813  if (value=='\n')
2814  line++;
2815  }
2816 
2817  if (value==EOF)
2818  error(line, "expected \"]\" in input file");
2819 }
2820 
2821 bool CHMM::comma_or_space(FILE* file)
2822 {
2823  int32_t value;
2824  while (((value=fgetc(file)) != EOF) && (value!=',') && (value!=';') && (value!=']')) //skip possible spaces and end if ',' or ';' occurs
2825  {
2826  if (value=='\n')
2827  line++;
2828  }
2829  if (value==']')
2830  {
2831  ungetc(value, file);
2832  SG_ERROR("found ']' instead of ';' or ','\n")
2833  return false ;
2834  } ;
2835 
2836  if (value==EOF)
2837  error(line, "expected \";\" or \",\" in input file");
2838 
2839  while (((value=fgetc(file)) != EOF) && (isspace(value))) //skip possible spaces
2840  {
2841  if (value=='\n')
2842  line++;
2843  }
2844  ungetc(value, file);
2845  return true ;
2846 }
2847 
2848 bool CHMM::get_numbuffer(FILE* file, char* buffer, int32_t length)
2849 {
2850  signed char value;
2851 
2852  while (((value=fgetc(file)) != EOF) &&
2853  !isdigit(value) && (value!='A')
2854  && (value!='C') && (value!='G') && (value!='T')
2855  && (value!='N') && (value!='n')
2856  && (value!='.') && (value!='-') && (value!='e') && (value!=']')) //skip possible spaces+crap
2857  {
2858  if (value=='\n')
2859  line++;
2860  }
2861  if (value==']')
2862  {
2863  ungetc(value,file) ;
2864  return false ;
2865  } ;
2866  if (value!=EOF)
2867  {
2868  int32_t i=0;
2869  switch (value)
2870  {
2871  case 'A':
2872  value='0' +CAlphabet::B_A;
2873  break;
2874  case 'C':
2875  value='0' +CAlphabet::B_C;
2876  break;
2877  case 'G':
2878  value='0' +CAlphabet::B_G;
2879  break;
2880  case 'T':
2881  value='0' +CAlphabet::B_T;
2882  break;
2883  };
2884 
2885  buffer[i++]=value;
2886 
2887  while (((value=fgetc(file)) != EOF) &&
2888  (isdigit(value) || (value=='.') || (value=='-') || (value=='e')
2889  || (value=='A') || (value=='C') || (value=='G')|| (value=='T')
2890  || (value=='N') || (value=='n')) && (i<length))
2891  {
2892  switch (value)
2893  {
2894  case 'A':
2895  value='0' +CAlphabet::B_A;
2896  break;
2897  case 'C':
2898  value='0' +CAlphabet::B_C;
2899  break;
2900  case 'G':
2901  value='0' +CAlphabet::B_G;
2902  break;
2903  case 'T':
2904  value='0' +CAlphabet::B_T;
2905  break;
2906  case '1': case '2': case'3': case '4': case'5':
2907  case '6': case '7': case'8': case '9': case '0': break ;
2908  case '.': case 'e': case '-': break ;
2909  default:
2910  SG_ERROR("found crap: %i %c (pos:%li)\n",i,value,ftell(file))
2911  };
2912  buffer[i++]=value;
2913  }
2914  ungetc(value, file);
2915  buffer[i]='\0';
2916 
2917  return (i<=length) && (i>0);
2918  }
2919  return false;
2920 }
2921 
2922 /*
2923  -format specs: model_file (model.hmm)
2924  % HMM - specification
2925  % N - number of states
2926  % M - number of observation_tokens
2927  % a is state_transition_matrix
2928  % size(a)= [N,N]
2929  %
2930  % b is observation_per_state_matrix
2931  % size(b)= [N,M]
2932  %
2933  % p is initial distribution
2934  % size(p)= [1, N]
2935 
2936  N=<int32_t>;
2937  M=<int32_t>;
2938 
2939  p=[<float64_t>,<float64_t>...<DOUBLE>];
2940  q=[<DOUBLE>,<DOUBLE>...<DOUBLE>];
2941 
2942  a=[ [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2943  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2944  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2945  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2946  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2947  ];
2948 
2949  b=[ [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2950  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2951  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2952  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2953  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
2954  ];
2955  */
2956 
2957 bool CHMM::load_model(FILE* file)
2958 {
2959  int32_t received_params=0; //a,b,p,N,M,O
2960 
2961  bool result=false;
2962  E_STATE state=INITIAL;
2963  char buffer[1024];
2964 
2965  line=1;
2966  int32_t i,j;
2967 
2968  if (file)
2969  {
2970  while (state!=END)
2971  {
2972  int32_t value=fgetc(file);
2973 
2974  if (value=='\n')
2975  line++;
2976  if (value==EOF)
2977  state=END;
2978 
2979  switch (state)
2980  {
2981  case INITIAL: // in the initial state only N,M initialisations and comments are allowed
2982  if (value=='N')
2983  {
2984  if (received_params & GOTN)
2985  error(line, "in model file: \"p double defined\"");
2986  else
2987  state=GET_N;
2988  }
2989  else if (value=='M')
2990  {
2991  if (received_params & GOTM)
2992  error(line, "in model file: \"p double defined\"");
2993  else
2994  state=GET_M;
2995  }
2996  else if (value=='%')
2997  {
2998  state=COMMENT;
2999  }
3000  case ARRAYs: // when n,m, order are known p,a,b arrays are allowed to be read
3001  if (value=='p')
3002  {
3003  if (received_params & GOTp)
3004  error(line, "in model file: \"p double defined\"");
3005  else
3006  state=GET_p;
3007  }
3008  if (value=='q')
3009  {
3010  if (received_params & GOTq)
3011  error(line, "in model file: \"q double defined\"");
3012  else
3013  state=GET_q;
3014  }
3015  else if (value=='a')
3016  {
3017  if (received_params & GOTa)
3018  error(line, "in model file: \"a double defined\"");
3019  else
3020  state=GET_a;
3021  }
3022  else if (value=='b')
3023  {
3024  if (received_params & GOTb)
3025  error(line, "in model file: \"b double defined\"");
3026  else
3027  state=GET_b;
3028  }
3029  else if (value=='%')
3030  {
3031  state=COMMENT;
3032  }
3033  break;
3034  case GET_N:
3035  if (value=='=')
3036  {
3037  if (get_numbuffer(file, buffer, 4)) //get num
3038  {
3039  this->N= atoi(buffer);
3040  received_params|=GOTN;
3041  state= (received_params == (GOTN | GOTM | GOTO)) ? ARRAYs : INITIAL;
3042  }
3043  else
3044  state=END; //end if error
3045  }
3046  break;
3047  case GET_M:
3048  if (value=='=')
3049  {
3050  if (get_numbuffer(file, buffer, 4)) //get num
3051  {
3052  this->M= atoi(buffer);
3053  received_params|=GOTM;
3054  state= (received_params == (GOTN | GOTM | GOTO)) ? ARRAYs : INITIAL;
3055  }
3056  else
3057  state=END; //end if error
3058  }
3059  break;
3060  case GET_a:
3061  if (value=='=')
3062  {
3063  float64_t f;
3064 
3065  transition_matrix_a=SG_MALLOC(float64_t, N*N);
3066  open_bracket(file);
3067  for (i=0; i<this->N; i++)
3068  {
3069  open_bracket(file);
3070 
3071  for (j=0; j<this->N ; j++)
3072  {
3073 
3074  if (fscanf( file, "%le", &f ) != 1)
3075  error(line, "float64_t expected");
3076  else
3077  set_a(i,j, f);
3078 
3079  if (j<this->N-1)
3080  comma_or_space(file);
3081  else
3082  close_bracket(file);
3083  }
3084 
3085  if (i<this->N-1)
3086  comma_or_space(file);
3087  else
3088  close_bracket(file);
3089  }
3090  received_params|=GOTa;
3091  }
3092  state= (received_params == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
3093  break;
3094  case GET_b:
3095  if (value=='=')
3096  {
3097  float64_t f;
3098 
3099  observation_matrix_b=SG_MALLOC(float64_t, N*M);
3100  open_bracket(file);
3101  for (i=0; i<this->N; i++)
3102  {
3103  open_bracket(file);
3104 
3105  for (j=0; j<this->M ; j++)
3106  {
3107 
3108  if (fscanf( file, "%le", &f ) != 1)
3109  error(line, "float64_t expected");
3110  else
3111  set_b(i,j, f);
3112 
3113  if (j<this->M-1)
3114  comma_or_space(file);
3115  else
3116  close_bracket(file);
3117  }
3118 
3119  if (i<this->N-1)
3120  comma_or_space(file);
3121  else
3122  close_bracket(file);
3123  }
3124  received_params|=GOTb;
3125  }
3126  state= ((received_params & (GOTa | GOTb | GOTp | GOTq)) == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
3127  break;
3128  case GET_p:
3129  if (value=='=')
3130  {
3131  float64_t f;
3132 
3133  initial_state_distribution_p=SG_MALLOC(float64_t, N);
3134  open_bracket(file);
3135  for (i=0; i<this->N ; i++)
3136  {
3137  if (fscanf( file, "%le", &f ) != 1)
3138  error(line, "float64_t expected");
3139  else
3140  set_p(i, f);
3141 
3142  if (i<this->N-1)
3143  comma_or_space(file);
3144  else
3145  close_bracket(file);
3146  }
3147  received_params|=GOTp;
3148  }
3149  state= (received_params == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
3150  break;
3151  case GET_q:
3152  if (value=='=')
3153  {
3154  float64_t f;
3155 
3156  end_state_distribution_q=SG_MALLOC(float64_t, N);
3157  open_bracket(file);
3158  for (i=0; i<this->N ; i++)
3159  {
3160  if (fscanf( file, "%le", &f ) != 1)
3161  error(line, "float64_t expected");
3162  else
3163  set_q(i, f);
3164 
3165  if (i<this->N-1)
3166  comma_or_space(file);
3167  else
3168  close_bracket(file);
3169  }
3170  received_params|=GOTq;
3171  }
3172  state= (received_params == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
3173  break;
3174  case COMMENT:
3175  if (value==EOF)
3176  state=END;
3177  else if (value=='\n')
3178  {
3179  line++;
3180  state=INITIAL;
3181  }
3182  break;
3183 
3184  default:
3185  break;
3186  }
3187  }
3188  result= (received_params== (GOTa | GOTb | GOTp | GOTq | GOTN | GOTM | GOTO));
3189  }
3190 
3191  SG_WARNING("not normalizing anymore, call normalize_hmm to make sure the hmm is valid!!\n")
3193  return result;
3194 }
3195 
3196 /*
3197  -format specs: train_file (train.trn)
3198  % HMM-TRAIN - specification
3199  % learn_a - elements in state_transition_matrix to be learned
3200  % learn_b - elements in oberservation_per_state_matrix to be learned
3201  % note: each line stands for
3202  % <state>, <observation(0)>, observation(1)...observation(NOW)>
3203  % learn_p - elements in initial distribution to be learned
3204  % learn_q - elements in the end-state distribution to be learned
3205  %
3206  % const_x - specifies initial values of elements
3207  % rest is assumed to be 0.0
3208  %
3209  % NOTE: IMPLICIT DEFINES:
3210  % #define A 0
3211  % #define C 1
3212  % #define G 2
3213  % #define T 3
3214  %
3215 
3216  learn_a=[ [<int32_t>,<int32_t>];
3217  [<int32_t>,<int32_t>];
3218  [<int32_t>,<int32_t>];
3219  ........
3220  [<int32_t>,<int32_t>];
3221  [-1,-1];
3222  ];
3223 
3224  learn_b=[ [<int32_t>,<int32_t>];
3225  [<int32_t>,<int32_t>];
3226  [<int32_t>,<int32_t>];
3227  ........
3228  [<int32_t>,<int32_t>];
3229  [-1,-1];
3230  ];
3231 
3232  learn_p= [ <int32_t>, ... , <int32_t>, -1 ];
3233  learn_q= [ <int32_t>, ... , <int32_t>, -1 ];
3234 
3235 
3236  const_a=[ [<int32_t>,<int32_t>,<DOUBLE>];
3237  [<int32_t>,<int32_t>,<DOUBLE>];
3238  [<int32_t>,<int32_t>,<DOUBLE>];
3239  ........
3240  [<int32_t>,<int32_t>,<DOUBLE>];
3241  [-1,-1,-1];
3242  ];
3243 
3244  const_b=[ [<int32_t>,<int32_t>,<DOUBLE>];
3245  [<int32_t>,<int32_t>,<DOUBLE>];
3246  [<int32_t>,<int32_t>,<DOUBLE];
3247  ........
3248  [<int32_t>,<int32_t>,<DOUBLE>];
3249  [-1,-1];
3250  ];
3251 
3252  const_p[]=[ [<int32_t>, <DOUBLE>], ... , [<int32_t>,<DOUBLE>], [-1,-1] ];
3253  const_q[]=[ [<int32_t>, <DOUBLE>], ... , [<int32_t>,<DOUBLE>], [-1,-1] ];
3254  */
3255 bool CHMM::load_definitions(FILE* file, bool verbose, bool _initialize)
3256 {
3257  if (model)
3258  delete model ;
3259  model=new Model();
3260 
3261  int32_t received_params=0x0000000; //a,b,p,q,N,M,O
3262  char buffer[1024];
3263 
3264  bool result=false;
3265  E_STATE state=INITIAL;
3266 
3267  { // do some useful initializations
3268  model->set_learn_a(0, -1);
3269  model->set_learn_a(1, -1);
3270  model->set_const_a(0, -1);
3271  model->set_const_a(1, -1);
3272  model->set_const_a_val(0, 1.0);
3273  model->set_learn_b(0, -1);
3274  model->set_const_b(0, -1);
3275  model->set_const_b_val(0, 1.0);
3276  model->set_learn_p(0, -1);
3277  model->set_learn_q(0, -1);
3278  model->set_const_p(0, -1);
3279  model->set_const_q(0, -1);
3280  } ;
3281 
3282  line=1;
3283 
3284  if (file)
3285  {
3286  while (state!=END)
3287  {
3288  int32_t value=fgetc(file);
3289 
3290  if (value=='\n')
3291  line++;
3292 
3293  if (value==EOF)
3294  state=END;
3295 
3296  switch (state)
3297  {
3298  case INITIAL:
3299  if (value=='l')
3300  {
3301  if (fgetc(file)=='e' && fgetc(file)=='a' && fgetc(file)=='r' && fgetc(file)=='n' && fgetc(file)=='_')
3302  {
3303  switch(fgetc(file))
3304  {
3305  case 'a':
3306  state=GET_learn_a;
3307  break;
3308  case 'b':
3309  state=GET_learn_b;
3310  break;
3311  case 'p':
3312  state=GET_learn_p;
3313  break;
3314  case 'q':
3315  state=GET_learn_q;
3316  break;
3317  default:
3318  error(line, "a,b,p or q expected in train definition file");
3319  };
3320  }
3321  }
3322  else if (value=='c')
3323  {
3324  if (fgetc(file)=='o' && fgetc(file)=='n' && fgetc(file)=='s'
3325  && fgetc(file)=='t' && fgetc(file)=='_')
3326  {
3327  switch(fgetc(file))
3328  {
3329  case 'a':
3330  state=GET_const_a;
3331  break;
3332  case 'b':
3333  state=GET_const_b;
3334  break;
3335  case 'p':
3336  state=GET_const_p;
3337  break;
3338  case 'q':
3339  state=GET_const_q;
3340  break;
3341  default:
3342  error(line, "a,b,p or q expected in train definition file");
3343  };
3344  }
3345  }
3346  else if (value=='%')
3347  {
3348  state=COMMENT;
3349  }
3350  else if (value==EOF)
3351  {
3352  state=END;
3353  }
3354  break;
3355  case GET_learn_a:
3356  if (value=='=')
3357  {
3358  open_bracket(file);
3359  bool finished=false;
3360  int32_t i=0;
3361 
3362  if (verbose)
3363  SG_DEBUG("\nlearn for transition matrix: ")
3364  while (!finished)
3365  {
3366  open_bracket(file);
3367 
3368  if (get_numbuffer(file, buffer, 4)) //get num
3369  {
3370  value=atoi(buffer);
3371  model->set_learn_a(i++, value);
3372 
3373  if (value<0)
3374  {
3375  finished=true;
3376  break;
3377  }
3378  if (value>=N)
3379  SG_ERROR("invalid value for learn_a(%i,0): %i\n",i/2,(int)value)
3380  }
3381  else
3382  break;
3383 
3384  comma_or_space(file);
3385 
3386  if (get_numbuffer(file, buffer, 4)) //get num
3387  {
3388  value=atoi(buffer);
3389  model->set_learn_a(i++, value);
3390 
3391  if (value<0)
3392  {
3393  finished=true;
3394  break;
3395  }
3396  if (value>=N)
3397  SG_ERROR("invalid value for learn_a(%i,1): %i\n",i/2-1,(int)value)
3398 
3399  }
3400  else
3401  break;
3402  close_bracket(file);
3403  }
3404  close_bracket(file);
3405  if (verbose)
3406  SG_DEBUG("%i Entries",(int)(i/2))
3407 
3408  if (finished)
3409  {
3410  received_params|=GOTlearn_a;
3411 
3412  state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
3413  }
3414  else
3415  state=END;
3416  }
3417  break;
3418  case GET_learn_b:
3419  if (value=='=')
3420  {
3421  open_bracket(file);
3422  bool finished=false;
3423  int32_t i=0;
3424 
3425  if (verbose)
3426  SG_DEBUG("\nlearn for emission matrix: ")
3427 
3428  while (!finished)
3429  {
3430  open_bracket(file);
3431 
3432  int32_t combine=0;
3433 
3434  for (int32_t j=0; j<2; j++)
3435  {
3436  if (get_numbuffer(file, buffer, 4)) //get num
3437  {
3438  value=atoi(buffer);
3439 
3440  if (j==0)
3441  {
3442  model->set_learn_b(i++, value);
3443 
3444  if (value<0)
3445  {
3446  finished=true;
3447  break;
3448  }
3449  if (value>=N)
3450  SG_ERROR("invalid value for learn_b(%i,0): %i\n",i/2,(int)value)
3451  }
3452  else
3453  combine=value;
3454  }
3455  else
3456  break;
3457 
3458  if (j<1)
3459  comma_or_space(file);
3460  else
3461  close_bracket(file);
3462  }
3463  model->set_learn_b(i++, combine);
3464  if (combine>=M)
3465 
3466  SG_ERROR("invalid value for learn_b(%i,1): %i\n",i/2-1,(int)value)
3467  }
3468  close_bracket(file);
3469  if (verbose)
3470  SG_DEBUG("%i Entries",(int)(i/2-1))
3471 
3472  if (finished)
3473  {
3474  received_params|=GOTlearn_b;
3475  state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
3476  }
3477  else
3478  state=END;
3479  }
3480  break;
3481  case GET_learn_p:
3482  if (value=='=')
3483  {
3484  open_bracket(file);
3485  bool finished=false;
3486  int32_t i=0;
3487 
3488  if (verbose)
3489  SG_DEBUG("\nlearn start states: ")
3490  while (!finished)
3491  {
3492  if (get_numbuffer(file, buffer, 4)) //get num
3493  {
3494  value=atoi(buffer);
3495 
3496  model->set_learn_p(i++, value);
3497 
3498  if (value<0)
3499  {
3500  finished=true;
3501  break;
3502  }
3503  if (value>=N)
3504  SG_ERROR("invalid value for learn_p(%i): %i\n",i-1,(int)value)
3505  }
3506  else
3507  break;
3508 
3509  comma_or_space(file);
3510  }
3511 
3512  close_bracket(file);
3513  if (verbose)
3514  SG_DEBUG("%i Entries",i-1)
3515 
3516  if (finished)
3517  {
3518  received_params|=GOTlearn_p;
3519  state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
3520  }
3521  else
3522  state=END;
3523  }
3524  break;
3525  case GET_learn_q:
3526  if (value=='=')
3527  {
3528  open_bracket(file);
3529  bool finished=false;
3530  int32_t i=0;
3531 
3532  if (verbose)
3533  SG_DEBUG("\nlearn terminal states: ")
3534  while (!finished)
3535  {
3536  if (get_numbuffer(file, buffer, 4)) //get num
3537  {
3538  value=atoi(buffer);
3539  model->set_learn_q(i++, value);
3540 
3541  if (value<0)
3542  {
3543  finished=true;
3544  break;
3545  }
3546  if (value>=N)
3547  SG_ERROR("invalid value for learn_q(%i): %i\n",i-1,(int)value)
3548  }
3549  else
3550  break;
3551 
3552  comma_or_space(file);
3553  }
3554 
3555  close_bracket(file);
3556  if (verbose)
3557  SG_DEBUG("%i Entries",i-1)
3558 
3559  if (finished)
3560  {
3561  received_params|=GOTlearn_q;
3562  state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
3563  }
3564  else
3565  state=END;
3566  }
3567  break;
3568  case GET_const_a:
3569  if (value=='=')
3570  {
3571  open_bracket(file);
3572  bool finished=false;
3573  int32_t i=0;
3574 
3575  if (verbose)
3576 #ifdef USE_HMMDEBUG
3577  SG_DEBUG("\nconst for transition matrix: \n")
3578 #else
3579  SG_DEBUG("\nconst for transition matrix: ")
3580 #endif
3581  while (!finished)
3582  {
3583  open_bracket(file);
3584 
3585  if (get_numbuffer(file, buffer, 4)) //get num
3586  {
3587  value=atoi(buffer);
3588  model->set_const_a(i++, value);
3589 
3590  if (value<0)
3591  {
3592  finished=true;
3593  model->set_const_a(i++, value);
3594  model->set_const_a_val((int32_t)i/2 - 1, value);
3595  break;
3596  }
3597  if (value>=N)
3598  SG_ERROR("invalid value for const_a(%i,0): %i\n",i/2,(int)value)
3599  }
3600  else
3601  break;
3602 
3603  comma_or_space(file);
3604 
3605  if (get_numbuffer(file, buffer, 4)) //get num
3606  {
3607  value=atoi(buffer);
3608  model->set_const_a(i++, value);
3609 
3610  if (value<0)
3611  {
3612  finished=true;
3613  model->set_const_a_val((int32_t)i/2 - 1, value);
3614  break;
3615  }
3616  if (value>=N)
3617  SG_ERROR("invalid value for const_a(%i,1): %i\n",i/2-1,(int)value)
3618  }
3619  else
3620  break;
3621 
3622  if (!comma_or_space(file))
3623  model->set_const_a_val((int32_t)i/2 - 1, 1.0);
3624  else
3625  if (get_numbuffer(file, buffer, 10)) //get num
3626  {
3627  float64_t dvalue=atof(buffer);
3628  model->set_const_a_val((int32_t)i/2 - 1, dvalue);
3629  if (dvalue<0)
3630  {
3631  finished=true;
3632  break;
3633  }
3634  if ((dvalue>1.0) || (dvalue<0.0))
3635  SG_ERROR("invalid value for const_a_val(%i): %e\n",(int)i/2-1,dvalue)
3636  }
3637  else
3638  model->set_const_a_val((int32_t)i/2 - 1, 1.0);
3639 
3640 #ifdef USE_HMMDEBUG
3641  if (verbose)
3642  SG_ERROR("const_a(%i,%i)=%e\n", model->get_const_a((int32_t)i/2-1,0),model->get_const_a((int32_t)i/2-1,1),model->get_const_a_val((int32_t)i/2-1))
3643 #endif
3644  close_bracket(file);
3645  }
3646  close_bracket(file);
3647  if (verbose)
3648  SG_DEBUG("%i Entries",(int)i/2-1)
3649 
3650  if (finished)
3651  {
3652  received_params|=GOTconst_a;
3653  state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
3654  }
3655  else
3656  state=END;
3657  }
3658  break;
3659 
3660  case GET_const_b:
3661  if (value=='=')
3662  {
3663  open_bracket(file);
3664  bool finished=false;
3665  int32_t i=0;
3666 
3667  if (verbose)
3668 #ifdef USE_HMMDEBUG
3669  SG_DEBUG("\nconst for emission matrix: \n")
3670 #else
3671  SG_DEBUG("\nconst for emission matrix: ")
3672 #endif
3673  while (!finished)
3674  {
3675  open_bracket(file);
3676  int32_t combine=0;
3677  for (int32_t j=0; j<3; j++)
3678  {
3679  if (get_numbuffer(file, buffer, 10)) //get num
3680  {
3681  if (j==0)
3682  {
3683  value=atoi(buffer);
3684 
3685  model->set_const_b(i++, value);
3686 
3687  if (value<0)
3688  {
3689  finished=true;
3690  //model->set_const_b_val((int32_t)(i-1)/2, value);
3691  break;
3692  }
3693  if (value>=N)
3694  SG_ERROR("invalid value for const_b(%i,0): %i\n",i/2-1,(int)value)
3695  }
3696  else if (j==2)
3697  {
3698  float64_t dvalue=atof(buffer);
3699  model->set_const_b_val((int32_t)(i-1)/2, dvalue);
3700  if (dvalue<0)
3701  {
3702  finished=true;
3703  break;
3704  } ;
3705  if ((dvalue>1.0) || (dvalue<0.0))
3706  SG_ERROR("invalid value for const_b_val(%i,1): %e\n",i/2-1,dvalue)
3707  }
3708  else
3709  {
3710  value=atoi(buffer);
3711  combine= value;
3712  } ;
3713  }
3714  else
3715  {
3716  if (j==2)
3717  model->set_const_b_val((int32_t)(i-1)/2, 1.0);
3718  break;
3719  } ;
3720  if (j<2)
3721  if ((!comma_or_space(file)) && (j==1))
3722  {
3723  model->set_const_b_val((int32_t)(i-1)/2, 1.0) ;
3724  break ;
3725  } ;
3726  }
3727  close_bracket(file);
3728  model->set_const_b(i++, combine);
3729  if (combine>=M)
3730  SG_ERROR("invalid value for const_b(%i,1): %i\n",i/2-1, combine)
3731 #ifdef USE_HMMDEBUG
3732  if (verbose && !finished)
3733  SG_ERROR("const_b(%i,%i)=%e\n", model->get_const_b((int32_t)i/2-1,0),model->get_const_b((int32_t)i/2-1,1),model->get_const_b_val((int32_t)i/2-1))
3734 #endif
3735  }
3736  close_bracket(file);
3737  if (verbose)
3738  SG_ERROR("%i Entries",(int)i/2-1)
3739 
3740  if (finished)
3741  {
3742  received_params|=GOTconst_b;
3743  state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
3744  }
3745  else
3746  state=END;
3747  }
3748  break;
3749  case GET_const_p:
3750  if (value=='=')
3751  {
3752  open_bracket(file);
3753  bool finished=false;
3754  int32_t i=0;
3755 
3756  if (verbose)
3757 #ifdef USE_HMMDEBUG
3758  SG_DEBUG("\nconst for start states: \n")
3759 #else
3760  SG_DEBUG("\nconst for start states: ")
3761 #endif
3762  while (!finished)
3763  {
3764  open_bracket(file);
3765 
3766  if (get_numbuffer(file, buffer, 4)) //get num
3767  {
3768  value=atoi(buffer);
3769  model->set_const_p(i, value);
3770 
3771  if (value<0)
3772  {
3773  finished=true;
3774  model->set_const_p_val(i++, value);
3775  break;
3776  }
3777  if (value>=N)
3778  SG_ERROR("invalid value for const_p(%i): %i\n",i,(int)value)
3779 
3780  }
3781  else
3782  break;
3783 
3784  if (!comma_or_space(file))
3785  model->set_const_p_val(i++, 1.0);
3786  else
3787  if (get_numbuffer(file, buffer, 10)) //get num
3788  {
3789  float64_t dvalue=atof(buffer);
3790  model->set_const_p_val(i++, dvalue);
3791  if (dvalue<0)
3792  {
3793  finished=true;
3794  break;
3795  }
3796  if ((dvalue>1) || (dvalue<0))
3797  SG_ERROR("invalid value for const_p_val(%i): %e\n",i,dvalue)
3798  }
3799  else
3800  model->set_const_p_val(i++, 1.0);
3801 
3802  close_bracket(file);
3803 
3804 #ifdef USE_HMMDEBUG
3805  if (verbose)
3806  SG_DEBUG("const_p(%i)=%e\n", model->get_const_p(i-1),model->get_const_p_val(i-1))
3807 #endif
3808  }
3809  if (verbose)
3810  SG_DEBUG("%i Entries",i-1)
3811 
3812  close_bracket(file);
3813 
3814  if (finished)
3815  {
3816  received_params|=GOTconst_p;
3817  state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
3818  }
3819  else
3820  state=END;
3821  }
3822  break;
3823  case GET_const_q:
3824  if (value=='=')
3825  {
3826  open_bracket(file);
3827  bool finished=false;
3828  if (verbose)
3829 #ifdef USE_HMMDEBUG
3830  SG_DEBUG("\nconst for terminal states: \n")
3831 #else
3832  SG_DEBUG("\nconst for terminal states: ")
3833 #endif
3834  int32_t i=0;
3835 
3836  while (!finished)
3837  {
3838  open_bracket(file) ;
3839  if (get_numbuffer(file, buffer, 4)) //get num
3840  {
3841  value=atoi(buffer);
3842  model->set_const_q(i, value);
3843  if (value<0)
3844  {
3845  finished=true;
3846  model->set_const_q_val(i++, value);
3847  break;
3848  }
3849  if (value>=N)
3850  SG_ERROR("invalid value for const_q(%i): %i\n",i,(int)value)
3851  }
3852  else
3853  break;
3854 
3855  if (!comma_or_space(file))
3856  model->set_const_q_val(i++, 1.0);
3857  else
3858  if (get_numbuffer(file, buffer, 10)) //get num
3859  {
3860  float64_t dvalue=atof(buffer);
3861  model->set_const_q_val(i++, dvalue);
3862  if (dvalue<0)
3863  {
3864  finished=true;
3865  break;
3866  }
3867  if ((dvalue>1) || (dvalue<0))
3868  SG_ERROR("invalid value for const_q_val(%i): %e\n",i,(double) dvalue)
3869  }
3870  else
3871  model->set_const_q_val(i++, 1.0);
3872 
3873  close_bracket(file);
3874 #ifdef USE_HMMDEBUG
3875  if (verbose)
3876  SG_DEBUG("const_q(%i)=%e\n", model->get_const_q(i-1),model->get_const_q_val(i-1))
3877 #endif
3878  }
3879  if (verbose)
3880  SG_DEBUG("%i Entries",i-1)
3881 
3882  close_bracket(file);
3883 
3884  if (finished)
3885  {
3886  received_params|=GOTconst_q;
3887  state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
3888  }
3889  else
3890  state=END;
3891  }
3892  break;
3893  case COMMENT:
3894  if (value==EOF)
3895  state=END;
3896  else if (value=='\n')
3897  state=INITIAL;
3898  break;
3899 
3900  default:
3901  break;
3902  }
3903  }
3904  }
3905 
3906  /*result=((received_params&(GOTlearn_a | GOTconst_a))!=0) ;
3907  result=result && ((received_params&(GOTlearn_b | GOTconst_b))!=0) ;
3908  result=result && ((received_params&(GOTlearn_p | GOTconst_p))!=0) ;
3909  result=result && ((received_params&(GOTlearn_q | GOTconst_q))!=0) ; */
3910  result=1 ;
3911  if (result)
3912  {
3913  model->sort_learn_a() ;
3914  model->sort_learn_b() ;
3915  if (_initialize)
3916  {
3917  init_model_defined(); ;
3918  convert_to_log();
3919  } ;
3920  }
3921  if (verbose)
3922  SG_DEBUG("\n")
3923  return result;
3924 }
3925 
3926 /*
3927  -format specs: model_file (model.hmm)
3928  % HMM - specification
3929  % N - number of states
3930  % M - number of observation_tokens
3931  % a is state_transition_matrix
3932  % size(a)= [N,N]
3933  %
3934  % b is observation_per_state_matrix
3935  % size(b)= [N,M]
3936  %
3937  % p is initial distribution
3938  % size(p)= [1, N]
3939 
3940  N=<int32_t>;
3941  M=<int32_t>;
3942 
3943  p=[<DOUBLE>,<DOUBLE>...<DOUBLE>];
3944 
3945  a=[ [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3946  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3947  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3948  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3949  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3950  ];
3951 
3952  b=[ [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3953  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3954  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3955  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3956  [<DOUBLE>,<DOUBLE>...<DOUBLE>];
3957  ];
3958  */
3959 
3960 bool CHMM::save_model(FILE* file)
3961 {
3962  bool result=false;
3963  int32_t i,j;
3964  const float32_t NAN_REPLACEMENT = (float32_t) CMath::ALMOST_NEG_INFTY ;
3965 
3966  if (file)
3967  {
3968  fprintf(file,"%s","% HMM - specification\n% N - number of states\n% M - number of observation_tokens\n% a is state_transition_matrix\n% size(a)= [N,N]\n%\n% b is observation_per_state_matrix\n% size(b)= [N,M]\n%\n% p is initial distribution\n% size(p)= [1, N]\n\n% q is distribution of end states\n% size(q)= [1, N]\n");
3969  fprintf(file,"N=%d;\n",N);
3970  fprintf(file,"M=%d;\n",M);
3971 
3972  fprintf(file,"p=[");
3973  for (i=0; i<N; i++)
3974  {
3975  if (i<N-1) {
3976  if (CMath::is_finite(get_p(i)))
3977  fprintf(file, "%e,", (double)get_p(i));
3978  else
3979  fprintf(file, "%f,", NAN_REPLACEMENT);
3980  }
3981  else {
3982  if (CMath::is_finite(get_p(i)))
3983  fprintf(file, "%e", (double)get_p(i));
3984  else
3985  fprintf(file, "%f", NAN_REPLACEMENT);
3986  }
3987  }
3988 
3989  fprintf(file,"];\n\nq=[");
3990  for (i=0; i<N; i++)
3991  {
3992  if (i<N-1) {
3993  if (CMath::is_finite(get_q(i)))
3994  fprintf(file, "%e,", (double)get_q(i));
3995  else
3996  fprintf(file, "%f,", NAN_REPLACEMENT);
3997  }
3998  else {
3999  if (CMath::is_finite(get_q(i)))
4000  fprintf(file, "%e", (double)get_q(i));
4001  else
4002  fprintf(file, "%f", NAN_REPLACEMENT);
4003  }
4004  }
4005  fprintf(file,"];\n\na=[");
4006 
4007  for (i=0; i<N; i++)
4008  {
4009  fprintf(file, "\t[");
4010 
4011  for (j=0; j<N; j++)
4012  {
4013  if (j<N-1) {
4014  if (CMath::is_finite(get_a(i,j)))
4015  fprintf(file, "%e,", (double)get_a(i,j));
4016  else
4017  fprintf(file, "%f,", NAN_REPLACEMENT);
4018  }
4019  else {
4020  if (CMath::is_finite(get_a(i,j)))
4021  fprintf(file, "%e];\n", (double)get_a(i,j));
4022  else
4023  fprintf(file, "%f];\n", NAN_REPLACEMENT);
4024  }
4025  }
4026  }
4027 
4028  fprintf(file," ];\n\nb=[");
4029 
4030  for (i=0; i<N; i++)
4031  {
4032  fprintf(file, "\t[");
4033 
4034  for (j=0; j<M; j++)
4035  {
4036  if (j<M-1) {
4037  if (CMath::is_finite(get_b(i,j)))
4038  fprintf(file, "%e,", (double)get_b(i,j));
4039  else
4040  fprintf(file, "%f,", NAN_REPLACEMENT);
4041  }
4042  else {
4043  if (CMath::is_finite(get_b(i,j)))
4044  fprintf(file, "%e];\n", (double)get_b(i,j));
4045  else
4046  fprintf(file, "%f];\n", NAN_REPLACEMENT);
4047  }
4048  }
4049  }
4050  result= (fprintf(file," ];\n") == 5);
4051  }
4052 
4053  return result;
4054 }
4055 
4056 T_STATES* CHMM::get_path(int32_t dim, float64_t& prob)
4057 {
4058  T_STATES* result = NULL;
4059 
4060  prob = best_path(dim);
4061  result = SG_MALLOC(T_STATES, p_observations->get_vector_length(dim));
4062 
4063  for (int32_t i=0; i<p_observations->get_vector_length(dim); i++)
4064  result[i]=PATH(dim)[i];
4065 
4066  return result;
4067 }
4068 
4069 bool CHMM::save_path(FILE* file)
4070 {
4071  bool result=false;
4072 
4073  if (file)
4074  {
4075  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
4076  {
4077  if (dim%100==0)
4078  SG_PRINT("%i..", dim)
4079  float64_t prob = best_path(dim);
4080  fprintf(file,"%i. path probability:%e\nstate sequence:\n", dim, prob);
4081  for (int32_t i=0; i<p_observations->get_vector_length(dim)-1; i++)
4082  fprintf(file,"%d ", PATH(dim)[i]);
4083  fprintf(file,"%d", PATH(dim)[p_observations->get_vector_length(dim)-1]);
4084  fprintf(file,"\n\n") ;
4085  }
4086  SG_DONE()
4087  result=true;
4088  }
4089 
4090  return result;
4091 }
4092 
4093 bool CHMM::save_likelihood_bin(FILE* file)
4094 {
4095  bool result=false;
4096 
4097  if (file)
4098  {
4099  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
4100  {
4101  float32_t prob= (float32_t) model_probability(dim);
4102  fwrite(&prob, sizeof(float32_t), 1, file);
4103  }
4104  result=true;
4105  }
4106 
4107  return result;
4108 }
4109 
4110 bool CHMM::save_likelihood(FILE* file)
4111 {
4112  bool result=false;
4113 
4114  if (file)
4115  {
4116  fprintf(file, "%% likelihood of model per observation\n%% P[O|model]=[ P[O|model]_1 P[O|model]_2 ... P[O|model]_dim ]\n%%\n");
4117 
4118  fprintf(file, "P=[");
4119  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
4120  fprintf(file, "%e ", (double) model_probability(dim));
4121 
4122  fprintf(file,"];");
4123  result=true;
4124  }
4125 
4126  return result;
4127 }
4128 
4129 #define FLOATWRITE(file, value) { float32_t rrr=float32_t(value); fwrite(&rrr, sizeof(float32_t), 1, file); num_floats++;}
4130 
4131 bool CHMM::save_model_bin(FILE* file)
4132 {
4133  int32_t i,j,q, num_floats=0 ;
4134  if (!model)
4135  {
4136  if (file)
4137  {
4138  // write id
4139  FLOATWRITE(file, (float32_t)CMath::INFTY);
4140  FLOATWRITE(file, (float32_t) 1);
4141 
4142  //derivates log(dp),log(dq)
4143  for (i=0; i<N; i++)
4144  FLOATWRITE(file, get_p(i));
4145  SG_INFO("wrote %i parameters for p\n",N)
4146 
4147  for (i=0; i<N; i++)
4148  FLOATWRITE(file, get_q(i)) ;
4149  SG_INFO("wrote %i parameters for q\n",N)
4150 
4151  //derivates log(da),log(db)
4152  for (i=0; i<N; i++)
4153  for (j=0; j<N; j++)
4154  FLOATWRITE(file, get_a(i,j));
4155  SG_INFO("wrote %i parameters for a\n",N*N)
4156 
4157  for (i=0; i<N; i++)
4158  for (j=0; j<M; j++)
4159  FLOATWRITE(file, get_b(i,j));
4160  SG_INFO("wrote %i parameters for b\n",N*M)
4161 
4162  // write id
4163  FLOATWRITE(file, (float32_t)CMath::INFTY);
4164  FLOATWRITE(file, (float32_t) 3);
4165 
4166  // write number of parameters
4167  FLOATWRITE(file, (float32_t) N);
4168  FLOATWRITE(file, (float32_t) N);
4169  FLOATWRITE(file, (float32_t) N*N);
4170  FLOATWRITE(file, (float32_t) N*M);
4171  FLOATWRITE(file, (float32_t) N);
4172  FLOATWRITE(file, (float32_t) M);
4173  } ;
4174  }
4175  else
4176  {
4177  if (file)
4178  {
4179  int32_t num_p, num_q, num_a, num_b ;
4180  // write id
4181  FLOATWRITE(file, (float32_t)CMath::INFTY);
4182  FLOATWRITE(file, (float32_t) 2);
4183 
4184  for (i=0; model->get_learn_p(i)>=0; i++)
4185  FLOATWRITE(file, get_p(model->get_learn_p(i)));
4186  num_p=i ;
4187  SG_INFO("wrote %i parameters for p\n",num_p)
4188 
4189  for (i=0; model->get_learn_q(i)>=0; i++)
4190  FLOATWRITE(file, get_q(model->get_learn_q(i)));
4191  num_q=i ;
4192  SG_INFO("wrote %i parameters for q\n",num_q)
4193 
4194  //derivates log(da),log(db)
4195  for (q=0; model->get_learn_a(q,1)>=0; q++)
4196  {
4197  i=model->get_learn_a(q,0) ;
4198  j=model->get_learn_a(q,1) ;
4199  FLOATWRITE(file, (float32_t)i);
4200  FLOATWRITE(file, (float32_t)j);
4201  FLOATWRITE(file, get_a(i,j));
4202  } ;
4203  num_a=q ;
4204  SG_INFO("wrote %i parameters for a\n",num_a)
4205 
4206  for (q=0; model->get_learn_b(q,0)>=0; q++)
4207  {
4208  i=model->get_learn_b(q,0) ;
4209  j=model->get_learn_b(q,1) ;
4210  FLOATWRITE(file, (float32_t)i);
4211  FLOATWRITE(file, (float32_t)j);
4212  FLOATWRITE(file, get_b(i,j));
4213  } ;
4214  num_b=q ;
4215  SG_INFO("wrote %i parameters for b\n",num_b)
4216 
4217  // write id
4218  FLOATWRITE(file, (float32_t)CMath::INFTY);
4219  FLOATWRITE(file, (float32_t) 3);
4220 
4221  // write number of parameters
4222  FLOATWRITE(file, (float32_t) num_p);
4223  FLOATWRITE(file, (float32_t) num_q);
4224  FLOATWRITE(file, (float32_t) num_a);
4225  FLOATWRITE(file, (float32_t) num_b);
4226  FLOATWRITE(file, (float32_t) N);
4227  FLOATWRITE(file, (float32_t) M);
4228  } ;
4229  } ;
4230  return true ;
4231 }
4232 
4233 bool CHMM::save_path_derivatives(FILE* logfile)
4234 {
4235  int32_t dim,i,j;
4236 
4237  if (logfile)
4238  {
4239  fprintf(logfile,"%% lambda denotes the model\n%% O denotes the observation sequence\n%% Q denotes the path\n%% \n%% calculating derivatives of P[O,Q|lambda]=p_{Q1}b_{Q1}(O_1}*a_{Q1}{Q2}b_{Q2}(O2)*...*q_{T-1}{T}b_{QT}(O_T}q_{Q_T} against p,q,a,b\n%%\n");
4240  fprintf(logfile,"%% dPr[...]=[ [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
4241  fprintf(logfile,"%% [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
4242  fprintf(logfile,"%% ............................. \n");
4243  fprintf(logfile,"%% [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_MM]\n");
4244  fprintf(logfile,"%% ];\n%%\n\ndPr(log()) = [\n");
4245  }
4246  else
4247  return false ;
4248 
4249  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
4250  {
4251  best_path(dim);
4252 
4253  fprintf(logfile, "[ ");
4254 
4255  //derivates dlogp,dlogq
4256  for (i=0; i<N; i++)
4257  fprintf(logfile,"%e, ", (double) path_derivative_p(i,dim) );
4258 
4259  for (i=0; i<N; i++)
4260  fprintf(logfile,"%e, ", (double) path_derivative_q(i,dim) );
4261 
4262  //derivates dloga,dlogb
4263  for (i=0; i<N; i++)
4264  for (j=0; j<N; j++)
4265  fprintf(logfile, "%e,", (double) path_derivative_a(i,j,dim) );
4266 
4267  for (i=0; i<N; i++)
4268  for (j=0; j<M; j++)
4269  fprintf(logfile, "%e,", (double) path_derivative_b(i,j,dim) );
4270 
4271  fseek(logfile,ftell(logfile)-1,SEEK_SET);
4272  fprintf(logfile, " ];\n");
4273  }
4274 
4275  fprintf(logfile, "];");
4276 
4277  return true ;
4278 }
4279 
4280 bool CHMM::save_path_derivatives_bin(FILE* logfile)
4281 {
4282  bool result=false;
4283  int32_t dim,i,j,q;
4284  float64_t prob=0 ;
4285  int32_t num_floats=0 ;
4286 
4287  float64_t sum_prob=0.0 ;
4288  if (!model)
4289  SG_WARNING("No definitions loaded -- writing derivatives of all weights\n")
4290  else
4291  SG_INFO("writing derivatives of changed weights only\n")
4292 
4293  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
4294  {
4295  if (dim%100==0)
4296  {
4297  SG_PRINT(".")
4298 
4299  } ;
4300 
4301  prob=best_path(dim);
4302  sum_prob+=prob ;
4303 
4304  if (!model)
4305  {
4306  if (logfile)
4307  {
4308  // write prob
4309  FLOATWRITE(logfile, prob);
4310 
4311  for (i=0; i<N; i++)
4312  FLOATWRITE(logfile, path_derivative_p(i,dim));
4313 
4314  for (i=0; i<N; i++)
4315  FLOATWRITE(logfile, path_derivative_q(i,dim));
4316 
4317  for (i=0; i<N; i++)
4318  for (j=0; j<N; j++)
4319  FLOATWRITE(logfile, path_derivative_a(i,j,dim));
4320 
4321  for (i=0; i<N; i++)
4322  for (j=0; j<M; j++)
4323  FLOATWRITE(logfile, path_derivative_b(i,j,dim));
4324 
4325  }
4326  }
4327  else
4328  {
4329  if (logfile)
4330  {
4331  // write prob
4332  FLOATWRITE(logfile, prob);
4333 
4334  for (i=0; model->get_learn_p(i)>=0; i++)
4335  FLOATWRITE(logfile, path_derivative_p(model->get_learn_p(i),dim));
4336 
4337  for (i=0; model->get_learn_q(i)>=0; i++)
4338  FLOATWRITE(logfile, path_derivative_q(model->get_learn_q(i),dim));
4339 
4340  for (q=0; model->get_learn_a(q,0)>=0; q++)
4341  {
4342  i=model->get_learn_a(q,0) ;
4343  j=model->get_learn_a(q,1) ;
4344  FLOATWRITE(logfile, path_derivative_a(i,j,dim));
4345  } ;
4346 
4347  for (q=0; model->get_learn_b(q,0)>=0; q++)
4348  {
4349  i=model->get_learn_b(q,0) ;
4350  j=model->get_learn_b(q,1) ;
4351  FLOATWRITE(logfile, path_derivative_b(i,j,dim));
4352  } ;
4353  }
4354  } ;
4355  }
4356  save_model_bin(logfile) ;
4357 
4358  result=true;
4359  SG_PRINT("\n")
4360  return result;
4361 }
4362 
4363 bool CHMM::save_model_derivatives_bin(FILE* file)
4364 {
4365  bool result=false;
4366  int32_t dim,i,j,q ;
4367  int32_t num_floats=0 ;
4368 
4369  if (!model)
4370  SG_WARNING("No definitions loaded -- writing derivatives of all weights\n")
4371  else
4372  SG_INFO("writing derivatives of changed weights only\n")
4373 
4374 #ifdef USE_HMMPARALLEL
4375  int32_t num_threads = parallel->get_num_threads();
4376  pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
4377  S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);
4378 
4379  if (p_observations->get_num_vectors()<num_threads)
4380  num_threads=p_observations->get_num_vectors();
4381 #endif
4382 
4383  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
4384  {
4385  if (dim%20==0)
4386  {
4387  SG_PRINT(".")
4388 
4389  } ;
4390 
4391 #ifdef USE_HMMPARALLEL
4392  if (dim%num_threads==0)
4393  {
4394  for (i=0; i<num_threads; i++)
4395  {
4396  if (dim+i<p_observations->get_num_vectors())
4397  {
4398  params[i].hmm=this ;
4399  params[i].dim=dim+i ;
4400  pthread_create(&threads[i], NULL, bw_dim_prefetch, (void*)&params[i]) ;
4401  }
4402  }
4403 
4404  for (i=0; i<num_threads; i++)
4405  {
4406  if (dim+i<p_observations->get_num_vectors())
4407  pthread_join(threads[i], NULL);
4408  }
4409  }
4410 #endif
4411 
4412  float64_t prob=model_probability(dim) ;
4413  if (!model)
4414  {
4415  if (file)
4416  {
4417  // write prob
4418  FLOATWRITE(file, prob);
4419 
4420  //derivates log(dp),log(dq)
4421  for (i=0; i<N; i++)
4422  FLOATWRITE(file, model_derivative_p(i,dim));
4423 
4424  for (i=0; i<N; i++)
4425  FLOATWRITE(file, model_derivative_q(i,dim));
4426 
4427  //derivates log(da),log(db)
4428  for (i=0; i<N; i++)
4429  for (j=0; j<N; j++)
4430  FLOATWRITE(file, model_derivative_a(i,j,dim));
4431 
4432  for (i=0; i<N; i++)
4433  for (j=0; j<M; j++)
4434  FLOATWRITE(file, model_derivative_b(i,j,dim));
4435 
4436  if (dim==0)
4437  SG_INFO("Number of parameters (including posterior prob.): %i\n", num_floats)
4438  } ;
4439  }
4440  else
4441  {
4442  if (file)
4443  {
4444  // write prob
4445  FLOATWRITE(file, prob);
4446 
4447  for (i=0; model->get_learn_p(i)>=0; i++)
4448  FLOATWRITE(file, model_derivative_p(model->get_learn_p(i),dim));
4449 
4450  for (i=0; model->get_learn_q(i)>=0; i++)
4451  FLOATWRITE(file, model_derivative_q(model->get_learn_q(i),dim));
4452 
4453  //derivates log(da),log(db)
4454  for (q=0; model->get_learn_a(q,1)>=0; q++)
4455  {
4456  i=model->get_learn_a(q,0) ;
4457  j=model->get_learn_a(q,1) ;
4458  FLOATWRITE(file, model_derivative_a(i,j,dim));
4459  } ;
4460 
4461  for (q=0; model->get_learn_b(q,0)>=0; q++)
4462  {
4463  i=model->get_learn_b(q,0) ;
4464  j=model->get_learn_b(q,1) ;
4465  FLOATWRITE(file, model_derivative_b(i,j,dim));
4466  } ;
4467  if (dim==0)
4468  SG_INFO("Number of parameters (including posterior prob.): %i\n", num_floats)
4469  } ;
4470  } ;
4471  }
4472  save_model_bin(file) ;
4473 
4474 #ifdef USE_HMMPARALLEL
4475  SG_FREE(threads);
4476  SG_FREE(params);
4477 #endif
4478 
4479  result=true;
4480  SG_PRINT("\n")
4481  return result;
4482 }
4483 
4484 bool CHMM::save_model_derivatives(FILE* file)
4485 {
4486  bool result=false;
4487  int32_t dim,i,j;
4488 
4489  if (file)
4490  {
4491 
4492  fprintf(file,"%% lambda denotes the model\n%% O denotes the observation sequence\n%% Q denotes the path\n%%\n%% calculating derivatives of P[O|lambda]=sum_{all Q}p_{Q1}b_{Q1}(O_1}*a_{Q1}{Q2}b_{Q2}(O2)*...*q_{T-1}{T}b_{QT}(O_T}q_{Q_T} against p,q,a,b\n%%\n");
4493  fprintf(file,"%% dPr[...]=[ [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
4494  fprintf(file,"%% [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
4495  fprintf(file,"%% ............................. \n");
4496  fprintf(file,"%% [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_MM]\n");
4497  fprintf(file,"%% ];\n%%\n\nlog(dPr) = [\n");
4498 
4499 
4500  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
4501  {
4502  fprintf(file, "[ ");
4503 
4504  //derivates log(dp),log(dq)
4505  for (i=0; i<N; i++)
4506  fprintf(file,"%e, ", (double) model_derivative_p(i, dim) ); //log (dp)
4507 
4508  for (i=0; i<N; i++)
4509  fprintf(file,"%e, ", (double) model_derivative_q(i, dim) ); //log (dq)
4510 
4511  //derivates log(da),log(db)
4512  for (i=0; i<N; i++)
4513  for (j=0; j<N; j++)
4514  fprintf(file, "%e,", (double) model_derivative_a(i,j,dim) );
4515 
4516  for (i=0; i<N; i++)
4517  for (j=0; j<M; j++)
4518  fprintf(file, "%e,", (double) model_derivative_b(i,j,dim) );
4519 
4520  fseek(file,ftell(file)-1,SEEK_SET);
4521  fprintf(file, " ];\n");
4522  }
4523 
4524 
4525  fprintf(file, "];");
4526 
4527  result=true;
4528  }
4529  return result;
4530 }
4531 
4532 bool CHMM::check_model_derivatives_combined()
4533 {
4534  // bool result=false;
4535  const float64_t delta=5e-4 ;
4536 
4537  int32_t i ;
4538  //derivates log(da)
4539  /* for (i=0; i<N; i++)
4540  {
4541  for (int32_t j=0; j<N; j++)
4542  {
4543  float64_t old_a=get_a(i,j) ;
4544 
4545  set_a(i,j, log(exp(old_a)-delta)) ;
4546  invalidate_model() ;
4547  float64_t prob_old=exp(model_probability(-1)*p_observations->get_num_vectors()) ;
4548 
4549  set_a(i,j, log(exp(old_a)+delta)) ;
4550  invalidate_model() ;
4551  float64_t prob_new=exp(model_probability(-1)*p_observations->get_num_vectors());
4552 
4553  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4554 
4555  set_a(i,j, old_a) ;
4556  invalidate_model() ;
4557 
4558  float64_t prod_prob=model_probability(-1)*p_observations->get_num_vectors() ;
4559 
4560  float64_t deriv_calc=0 ;
4561  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
4562  deriv_calc+=exp(model_derivative_a(i, j, dim)+
4563  prod_prob-model_probability(dim)) ;
4564 
4565  SG_DEBUG("da(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
4566  } ;
4567  } ;*/
4568  //derivates log(db)
4569  i=0;//for (i=0; i<N; i++)
4570  {
4571  for (int32_t j=0; j<M; j++)
4572  {
4573  float64_t old_b=get_b(i,j) ;
4574 
4575  set_b(i,j, log(exp(old_b)-delta)) ;
4576  invalidate_model() ;
4577  float64_t prob_old=(model_probability(-1)*p_observations->get_num_vectors()) ;
4578 
4579  set_b(i,j, log(exp(old_b)+delta)) ;
4580  invalidate_model() ;
4581  float64_t prob_new=(model_probability(-1)*p_observations->get_num_vectors());
4582 
4583  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4584 
4585  set_b(i,j, old_b) ;
4586  invalidate_model() ;
4587 
4588  float64_t deriv_calc=0 ;
4589  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
4590  {
4591  deriv_calc+=exp(model_derivative_b(i, j, dim)-model_probability(dim)) ;
4592  if (j==1)
4593  SG_INFO("deriv_calc[%i]=%e\n",dim,deriv_calc)
4594  } ;
4595 
4596  SG_ERROR("b(%i,%i)=%e db(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j,exp(old_b),i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
4597  } ;
4598  } ;
4599  return true ;
4600 }
4601 
4602 bool CHMM::check_model_derivatives()
4603 {
4604  bool result=false;
4605  const float64_t delta=3e-4 ;
4606 
4607  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
4608  {
4609  int32_t i ;
4610  //derivates log(dp),log(dq)
4611  for (i=0; i<N; i++)
4612  {
4613  for (int32_t j=0; j<N; j++)
4614  {
4615  float64_t old_a=get_a(i,j) ;
4616 
4617  set_a(i,j, log(exp(old_a)-delta)) ;
4618  invalidate_model() ;
4619  float64_t prob_old=exp(model_probability(dim)) ;
4620 
4621  set_a(i,j, log(exp(old_a)+delta)) ;
4622  invalidate_model() ;
4623  float64_t prob_new=exp(model_probability(dim));
4624 
4625  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4626 
4627  set_a(i,j, old_a) ;
4628  invalidate_model() ;
4629  float64_t deriv_calc=exp(model_derivative_a(i, j, dim)) ;
4630 
4631  SG_DEBUG("da(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
4632  invalidate_model() ;
4633  } ;
4634  } ;
4635  for (i=0; i<N; i++)
4636  {
4637  for (int32_t j=0; j<M; j++)
4638  {
4639  float64_t old_b=get_b(i,j) ;
4640 
4641  set_b(i,j, log(exp(old_b)-delta)) ;
4642  invalidate_model() ;
4643  float64_t prob_old=exp(model_probability(dim)) ;
4644 
4645  set_b(i,j, log(exp(old_b)+delta)) ;
4646  invalidate_model() ;
4647  float64_t prob_new=exp(model_probability(dim));
4648 
4649  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4650 
4651  set_b(i,j, old_b) ;
4652  invalidate_model() ;
4653  float64_t deriv_calc=exp(model_derivative_b(i, j, dim));
4654 
4655  SG_DEBUG("db(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/(deriv_calc))
4656  } ;
4657  } ;
4658 
4659 #ifdef TEST
4660  for (i=0; i<N; i++)
4661  {
4662  float64_t old_p=get_p(i) ;
4663 
4664  set_p(i, log(exp(old_p)-delta)) ;
4665  invalidate_model() ;
4666  float64_t prob_old=exp(model_probability(dim)) ;
4667 
4668  set_p(i, log(exp(old_p)+delta)) ;
4669  invalidate_model() ;
4670  float64_t prob_new=exp(model_probability(dim));
4671  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4672 
4673  set_p(i, old_p) ;
4674  invalidate_model() ;
4675  float64_t deriv_calc=exp(model_derivative_p(i, dim));
4676 
4677  //if (fabs(deriv_calc_old-deriv)>1e-4)
4678  SG_DEBUG("dp(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
4679  } ;
4680  for (i=0; i<N; i++)
4681  {
4682  float64_t old_q=get_q(i) ;
4683 
4684  set_q(i, log(exp(old_q)-delta)) ;
4685  invalidate_model() ;
4686  float64_t prob_old=exp(model_probability(dim)) ;
4687 
4688  set_q(i, log(exp(old_q)+delta)) ;
4689  invalidate_model() ;
4690  float64_t prob_new=exp(model_probability(dim));
4691 
4692  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4693 
4694  set_q(i, old_q) ;
4695  invalidate_model() ;
4696  float64_t deriv_calc=exp(model_derivative_q(i, dim));
4697 
4698  //if (fabs(deriv_calc_old-deriv)>1e-4)
4699  SG_DEBUG("dq(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
4700  } ;
4701 #endif
4702  }
4703  return result;
4704 }
4705 
4706 #ifdef USE_HMMDEBUG
4707 bool CHMM::check_path_derivatives()
4708 {
4709  bool result=false;
4710  const float64_t delta=1e-4 ;
4711 
4712  for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
4713  {
4714  int32_t i ;
4715  //derivates log(dp),log(dq)
4716  for (i=0; i<N; i++)
4717  {
4718  for (int32_t j=0; j<N; j++)
4719  {
4720  float64_t old_a=get_a(i,j) ;
4721 
4722  set_a(i,j, log(exp(old_a)-delta)) ;
4723  invalidate_model() ;
4724  float64_t prob_old=best_path(dim) ;
4725 
4726  set_a(i,j, log(exp(old_a)+delta)) ;
4727  invalidate_model() ;
4728  float64_t prob_new=best_path(dim);
4729 
4730  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4731 
4732  set_a(i,j, old_a) ;
4733  invalidate_model() ;
4734  float64_t deriv_calc=path_derivative_a(i, j, dim) ;
4735 
4736  SG_DEBUG("da(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
4737  } ;
4738  } ;
4739  for (i=0; i<N; i++)
4740  {
4741  for (int32_t j=0; j<M; j++)
4742  {
4743  float64_t old_b=get_b(i,j) ;
4744 
4745  set_b(i,j, log(exp(old_b)-delta)) ;
4746  invalidate_model() ;
4747  float64_t prob_old=best_path(dim) ;
4748 
4749  set_b(i,j, log(exp(old_b)+delta)) ;
4750  invalidate_model() ;
4751  float64_t prob_new=best_path(dim);
4752 
4753  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4754 
4755  set_b(i,j, old_b) ;
4756  invalidate_model() ;
4757  float64_t deriv_calc=path_derivative_b(i, j, dim);
4758 
4759  SG_DEBUG("db(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/(deriv_calc))
4760  } ;
4761  } ;
4762 
4763  for (i=0; i<N; i++)
4764  {
4765  float64_t old_p=get_p(i) ;
4766 
4767  set_p(i, log(exp(old_p)-delta)) ;
4768  invalidate_model() ;
4769  float64_t prob_old=best_path(dim) ;
4770 
4771  set_p(i, log(exp(old_p)+delta)) ;
4772  invalidate_model() ;
4773  float64_t prob_new=best_path(dim);
4774  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4775 
4776  set_p(i, old_p) ;
4777  invalidate_model() ;
4778  float64_t deriv_calc=path_derivative_p(i, dim);
4779 
4780  //if (fabs(deriv_calc_old-deriv)>1e-4)
4781  SG_DEBUG("dp(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
4782  } ;
4783  for (i=0; i<N; i++)
4784  {
4785  float64_t old_q=get_q(i) ;
4786 
4787  set_q(i, log(exp(old_q)-delta)) ;
4788  invalidate_model() ;
4789  float64_t prob_old=best_path(dim) ;
4790 
4791  set_q(i, log(exp(old_q)+delta)) ;
4792  invalidate_model() ;
4793  float64_t prob_new=best_path(dim);
4794 
4795  float64_t deriv = (prob_new-prob_old)/(2*delta) ;
4796 
4797  set_q(i, old_q) ;
4798  invalidate_model() ;
4799  float64_t deriv_calc=path_derivative_q(i, dim);
4800 
4801  //if (fabs(deriv_calc_old-deriv)>1e-4)
4802  SG_DEBUG("dq(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
4803  } ;
4804  }
4805  return result;
4806 }
4807 #endif // USE_HMMDEBUG
4808 
4809 //normalize model (sum to one constraint)
4810 void CHMM::normalize(bool keep_dead_states)
4811 {
4812  int32_t i,j;
4813  const float64_t INF=-1e10;
4814  float64_t sum_p =INF;
4815 
4816  for (i=0; i<N; i++)
4817  {
4818  sum_p=CMath::logarithmic_sum(sum_p, get_p(i));
4819 
4820  float64_t sum_b =INF;
4821  float64_t sum_a =get_q(i);
4822 
4823  for (j=0; j<N; j++)
4824  sum_a=CMath::logarithmic_sum(sum_a, get_a(i,j));
4825 
4826  if (sum_a>CMath::ALMOST_NEG_INFTY/N || (!keep_dead_states) )
4827  {
4828  for (j=0; j<N; j++)
4829  set_a(i,j, get_a(i,j)-sum_a);
4830  set_q(i, get_q(i)-sum_a);
4831  }
4832 
4833  for (j=0; j<M; j++)
4834  sum_b=CMath::logarithmic_sum(sum_b, get_b(i,j));
4835  for (j=0; j<M; j++)
4836  set_b(i,j, get_b(i,j)-sum_b);
4837  }
4838 
4839  for (i=0; i<N; i++)
4840  set_p(i, get_p(i)-sum_p);
4841 
4842  invalidate_model();
4843 }
4844 
4845 bool CHMM::append_model(CHMM* app_model)
4846 {
4847  bool result=false;
4848  const int32_t num_states=app_model->get_N();
4849  int32_t i,j;
4850 
4851  SG_DEBUG("cur N:%d M:%d\n", N, M)
4852  SG_DEBUG("old N:%d M:%d\n", app_model->get_N(), app_model->get_M())
4853  if (app_model->get_M() == get_M())
4854  {
4855  float64_t* n_p=SG_MALLOC(float64_t, N+num_states);
4856  float64_t* n_q=SG_MALLOC(float64_t, N+num_states);
4857  float64_t* n_a=SG_MALLOC(float64_t, (N+num_states)*(N+num_states));
4858  //SG_PRINT("size n_b: %d\n", (N+num_states)*M)
4859  float64_t* n_b=SG_MALLOC(float64_t, (N+num_states)*M);
4860 
4861  //clear n_x
4862  for (i=0; i<N+num_states; i++)
4863  {
4864  n_p[i]=-CMath::INFTY;
4865  n_q[i]=-CMath::INFTY;
4866 
4867  for (j=0; j<N+num_states; j++)
4868  n_a[(N+num_states)*i+j]=-CMath::INFTY;
4869 
4870  for (j=0; j<M; j++)
4871  n_b[M*i+j]=-CMath::INFTY;
4872  }
4873 
4874  //copy models first
4875  // warning pay attention to the ordering of
4876  // transition_matrix_a, observation_matrix_b !!!
4877 
4878  // cur_model
4879  for (i=0; i<N; i++)
4880  {
4881  n_p[i]=get_p(i);
4882 
4883  for (j=0; j<N; j++)
4884  n_a[(N+num_states)*j+i]=get_a(i,j);
4885 
4886  for (j=0; j<M; j++)
4887  {
4888  n_b[M*i+j]=get_b(i,j);
4889  }
4890  }
4891 
4892  // append_model
4893  for (i=0; i<app_model->get_N(); i++)
4894  {
4895  n_q[i+N]=app_model->get_q(i);
4896 
4897  for (j=0; j<app_model->get_N(); j++)
4898  n_a[(N+num_states)*(j+N)+(i+N)]=app_model->get_a(i,j);
4899  for (j=0; j<app_model->get_M(); j++)
4900  n_b[M*(i+N)+j]=app_model->get_b(i,j);
4901  }
4902 
4903 
4904  // transition to the two and back
4905  for (i=0; i<N; i++)
4906  {
4907  for (j=N; j<N+num_states; j++)
4908  n_a[(N+num_states)*j + i]=CMath::logarithmic_sum(get_q(i)+app_model->get_p(j-N), n_a[(N+num_states)*j + i]);
4909  }
4910 
4911  free_state_dependend_arrays();
4912  N+=num_states;
4913 
4914  alloc_state_dependend_arrays();
4915 
4916  //delete + adjust pointers
4917  SG_FREE(initial_state_distribution_p);
4918  SG_FREE(end_state_distribution_q);
4919  SG_FREE(transition_matrix_a);
4920  SG_FREE(observation_matrix_b);
4921 
4922  transition_matrix_a=n_a;
4923  observation_matrix_b=n_b;
4924  initial_state_distribution_p=n_p;
4925  end_state_distribution_q=n_q;
4926 
4927  SG_WARNING("not normalizing anymore, call normalize_hmm to make sure the hmm is valid!!\n")
4929  invalidate_model();
4930  }
4931  else
4932  SG_ERROR("number of observations is different for append model, doing nothing!\n")
4933 
4934  return result;
4935 }
4936 
4937 bool CHMM::append_model(CHMM* app_model, float64_t* cur_out, float64_t* app_out)
4938 {
4939  bool result=false;
4940  const int32_t num_states=app_model->get_N()+2;
4941  int32_t i,j;
4942 
4943  if (app_model->get_M() == get_M())
4944  {
4945  float64_t* n_p=SG_MALLOC(float64_t, N+num_states);
4946  float64_t* n_q=SG_MALLOC(float64_t, N+num_states);
4947  float64_t* n_a=SG_MALLOC(float64_t, (N+num_states)*(N+num_states));
4948  //SG_PRINT("size n_b: %d\n", (N+num_states)*M)
4949  float64_t* n_b=SG_MALLOC(float64_t, (N+num_states)*M);
4950 
4951  //clear n_x
4952  for (i=0; i<N+num_states; i++)
4953  {
4954  n_p[i]=-CMath::INFTY;
4955  n_q[i]=-CMath::INFTY;
4956 
4957  for (j=0; j<N+num_states; j++)
4958  n_a[(N+num_states)*j+i]=-CMath::INFTY;
4959 
4960  for (j=0; j<M; j++)
4961  n_b[M*i+j]=-CMath::INFTY;
4962  }
4963 
4964  //copy models first
4965  // warning pay attention to the ordering of
4966  // transition_matrix_a, observation_matrix_b !!!
4967 
4968  // cur_model
4969  for (i=0; i<N; i++)
4970  {
4971  n_p[i]=get_p(i);
4972 
4973  for (j=0; j<N; j++)
4974  n_a[(N+num_states)*j+i]=get_a(i,j);
4975 
4976  for (j=0; j<M; j++)
4977  {
4978  n_b[M*i+j]=get_b(i,j);
4979  }
4980  }
4981 
4982  // append_model
4983  for (i=0; i<app_model->get_N(); i++)
4984  {
4985  n_q[i+N+2]=app_model->get_q(i);
4986 
4987  for (j=0; j<app_model->get_N(); j++)
4988  n_a[(N+num_states)*(j+N+2)+(i+N+2)]=app_model->get_a(i,j);
4989  for (j=0; j<app_model->get_M(); j++)
4990  n_b[M*(i+N+2)+j]=app_model->get_b(i,j);
4991  }
4992 
4993  //initialize the two special states
4994 
4995  // output
4996  for (i=0; i<M; i++)
4997  {
4998  n_b[M*N+i]=cur_out[i];
4999  n_b[M*(N+1)+i]=app_out[i];
5000  }
5001 
5002  // transition to the two and back
5003  for (i=0; i<N+num_states; i++)
5004  {
5005  // the first state is only connected to the second
5006  if (i==N+1)
5007  n_a[(N+num_states)*i + N]=0;
5008 
5009  // only states of the cur_model can reach the
5010  // first state
5011  if (i<N)
5012  n_a[(N+num_states)*N+i]=get_q(i);
5013 
5014  // the second state is only connected to states of
5015  // the append_model (with probab app->p(i))
5016  if (i>=N+2)
5017  n_a[(N+num_states)*i+(N+1)]=app_model->get_p(i-(N+2));
5018  }
5019 
5020  free_state_dependend_arrays();
5021  N+=num_states;
5022 
5023  alloc_state_dependend_arrays();
5024 
5025  //delete + adjust pointers
5026  SG_FREE(initial_state_distribution_p);
5027  SG_FREE(end_state_distribution_q);
5028  SG_FREE(transition_matrix_a);
5029  SG_FREE(observation_matrix_b);
5030 
5031  transition_matrix_a=n_a;
5032  observation_matrix_b=n_b;
5033  initial_state_distribution_p=n_p;
5034  end_state_distribution_q=n_q;
5035 
5036  SG_WARNING("not normalizing anymore, call normalize_hmm to make sure the hmm is valid!!\n")
5038  invalidate_model();
5039  }
5040 
5041  return result;
5042 }
5043 
5044 
5045 void CHMM::add_states(int32_t num_states, float64_t default_value)
5046 {
5047  int32_t i,j;
5048  const float64_t MIN_RAND=1e-2; //this is the range of the random values for the new variables
5049  const float64_t MAX_RAND=2e-1;
5050 
5051  float64_t* n_p=SG_MALLOC(float64_t, N+num_states);
5052  float64_t* n_q=SG_MALLOC(float64_t, N+num_states);
5053  float64_t* n_a=SG_MALLOC(float64_t, (N+num_states)*(N+num_states));
5054  //SG_PRINT("size n_b: %d\n", (N+num_states)*M)
5055  float64_t* n_b=SG_MALLOC(float64_t, (N+num_states)*M);
5056 
5057  // warning pay attention to the ordering of
5058  // transition_matrix_a, observation_matrix_b !!!
5059  for (i=0; i<N; i++)
5060  {
5061  n_p[i]=get_p(i);
5062  n_q[i]=get_q(i);
5063 
5064  for (j=0; j<N; j++)
5065  n_a[(N+num_states)*j+i]=get_a(i,j);
5066 
5067  for (j=0; j<M; j++)
5068  n_b[M*i+j]=get_b(i,j);
5069  }
5070 
5071  for (i=N; i<N+num_states; i++)
5072  {
5073  n_p[i]=VAL_MACRO;
5074  n_q[i]=VAL_MACRO;
5075 
5076  for (j=0; j<N; j++)
5077  n_a[(N+num_states)*i+j]=VAL_MACRO;
5078 
5079  for (j=0; j<N+num_states; j++)
5080  n_a[(N+num_states)*j+i]=VAL_MACRO;
5081 
5082  for (j=0; j<M; j++)
5083  n_b[M*i+j]=VAL_MACRO;
5084  }
5085  free_state_dependend_arrays();
5086  N+=num_states;
5087 
5088  alloc_state_dependend_arrays();
5089 
5090  //delete + adjust pointers
5091  SG_FREE(initial_state_distribution_p);
5092  SG_FREE(end_state_distribution_q);
5093  SG_FREE(transition_matrix_a);
5094  SG_FREE(observation_matrix_b);
5095 
5096  transition_matrix_a=n_a;
5097  observation_matrix_b=n_b;
5098  initial_state_distribution_p=n_p;
5099  end_state_distribution_q=n_q;
5100 
5101  invalidate_model();
5102  normalize();
5103 }
5104 
5105 void CHMM::chop(float64_t value)
5106 {
5107  for (int32_t i=0; i<N; i++)
5108  {
5109  int32_t j;
5110 
5111  if (exp(get_p(i)) < value)
5112  set_p(i, CMath::ALMOST_NEG_INFTY);
5113 
5114  if (exp(get_q(i)) < value)
5115  set_q(i, CMath::ALMOST_NEG_INFTY);
5116 
5117  for (j=0; j<N; j++)
5118  {
5119  if (exp(get_a(i,j)) < value)
5120  set_a(i,j, CMath::ALMOST_NEG_INFTY);
5121  }
5122 
5123  for (j=0; j<M; j++)
5124  {
5125  if (exp(get_b(i,j)) < value)
5126  set_b(i,j, CMath::ALMOST_NEG_INFTY);
5127  }
5128  }
5129  normalize();
5130  invalidate_model();
5131 }
5132 
5133 bool CHMM::linear_train(bool right_align)
5134 {
5135  if (p_observations)
5136  {
5137  int32_t histsize=(get_M()*get_N());
5138  int32_t* hist=SG_MALLOC(int32_t, histsize);
5139  int32_t* startendhist=SG_MALLOC(int32_t, get_N());
5140  int32_t i,dim;
5141 
5142  ASSERT(p_observations->get_max_vector_length()<=get_N())
5143 
5144  for (i=0; i<histsize; i++)
5145  hist[i]=0;
5146 
5147  for (i=0; i<get_N(); i++)
5148  startendhist[i]=0;
5149 
5150  if (right_align)
5151  {
5152  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
5153  {
5154  int32_t len=0;
5155  bool free_vec;
5156  uint16_t* obs=p_observations->get_feature_vector(dim, len, free_vec);
5157 
5158  ASSERT(len<=get_N())
5159  startendhist[(get_N()-len)]++;
5160 
5161  for (i=0;i<len;i++)
5162  hist[(get_N()-len+i)*get_M() + *obs++]++;
5163 
5164  p_observations->free_feature_vector(obs, dim, free_vec);
5165  }
5166 
5167  set_q(get_N()-1, 1);
5168  for (i=0; i<get_N()-1; i++)
5169  set_q(i, 0);
5170 
5171  for (i=0; i<get_N(); i++)
5172  set_p(i, startendhist[i]+PSEUDO);
5173 
5174  for (i=0;i<get_N();i++)
5175  {
5176  for (int32_t j=0; j<get_N(); j++)
5177  {
5178  if (i==j-1)
5179  set_a(i,j, 1);
5180  else
5181  set_a(i,j, 0);
5182  }
5183  }
5184  }
5185  else
5186  {
5187  for (dim=0; dim<p_observations->get_num_vectors(); dim++)
5188  {
5189  int32_t len=0;
5190  bool free_vec;
5191  uint16_t* obs=p_observations->get_feature_vector(dim, len, free_vec);
5192 
5193  ASSERT(len<=get_N())
5194  for (i=0;i<len;i++)
5195  hist[i*get_M() + *obs++]++;
5196 
5197  startendhist[len-1]++;
5198 
5199  p_observations->free_feature_vector(obs, dim, free_vec);
5200  }
5201 
5202  set_p(0, 1);
5203  for (i=1; i<get_N(); i++)
5204  set_p(i, 0);
5205 
5206  for (i=0; i<get_N(); i++)
5207  set_q(i, startendhist[i]+PSEUDO);
5208 
5209  int32_t total=p_observations->get_num_vectors();
5210 
5211  for (i=0;i<get_N();i++)
5212  {
5213  total-= startendhist[i] ;
5214 
5215  for (int32_t j=0; j<get_N(); j++)
5216  {
5217  if (i==j-1)
5218  set_a(i,j, total+PSEUDO);
5219  else
5220  set_a(i,j, 0);
5221  }
5222  }
5223  ASSERT(total==0)
5224  }
5225 
5226  for (i=0;i<get_N();i++)
5227  {
5228  for (int32_t j=0; j<get_M(); j++)
5229  {
5230  float64_t sum=0;
5231  int32_t offs=i*get_M()+ p_observations->get_masked_symbols((uint16_t) j, (uint8_t) 254);
5232 
5233  for (int32_t k=0; k<p_observations->get_original_num_symbols(); k++)
5234  sum+=hist[offs+k];
5235 
5236  set_b(i,j, (PSEUDO+hist[i*get_M()+j])/(sum+PSEUDO*p_observations->get_original_num_symbols()));
5237  }
5238  }
5239 
5240  SG_FREE(hist);
5241  SG_FREE(startendhist);
5242  convert_to_log();
5243  invalidate_model();
5244  return true;
5245  }
5246  else
5247  return false;
5248 }
5249 
5250 void CHMM::set_observation_nocache(CStringFeatures<uint16_t>* obs)
5251 {
5252  ASSERT(obs)
5253  p_observations=obs;
5254  SG_REF(obs);
5255 
5256  if (obs)
5257  if (obs->get_num_symbols() > M)
5258  SG_ERROR("number of symbols in observation (%ld) larger than M (%d)\n", (long) obs->get_num_symbols(), M)
5259 
5260  if (!reused_caches)
5261  {
5262 #ifdef USE_HMMPARALLEL_STRUCTURES
5263  for (int32_t i=0; i<parallel->get_num_threads(); i++)
5264  {
5265  SG_FREE(alpha_cache[i].table);
5266  SG_FREE(beta_cache[i].table);
5267  SG_FREE(states_per_observation_psi[i]);
5268  SG_FREE(path[i]);
5269 
5270  alpha_cache[i].table=NULL;
5271  beta_cache[i].table=NULL;
5272  states_per_observation_psi[i]=NULL;
5273  path[i]=NULL;
5274  } ;
5275 #else
5276  SG_FREE(alpha_cache.table);
5277  SG_FREE(beta_cache.table);
5278  SG_FREE(states_per_observation_psi);
5279  SG_FREE(path);
5280 
5281  alpha_cache.table=NULL;
5282  beta_cache.table=NULL;
5283  states_per_observation_psi=NULL;
5284  path=NULL;
5285 
5286 #endif //USE_HMMPARALLEL_STRUCTURES
5287  }
5288 
5289  invalidate_model();
5290 }
5291 
5292 void CHMM::set_observations(CStringFeatures<uint16_t>* obs, CHMM* lambda)
5293 {
5294  ASSERT(obs)
5295  SG_REF(obs);
5296  p_observations=obs;
5297 
5298  /* from Distribution, necessary for calls to base class methods, like
5299  * get_log_likelihood_sample():
5300  */
5301  SG_REF(obs);
5302  features=obs;
5303 
5304  SG_DEBUG("num symbols alphabet: %ld\n", obs->get_alphabet()->get_num_symbols())
5305  SG_DEBUG("num symbols: %ld\n", obs->get_num_symbols())
5306  SG_DEBUG("M: %d\n", M)
5307 
5308  if (obs)
5309  {
5310  if (obs->get_num_symbols() > M)
5311  {
5312  SG_ERROR("number of symbols in observation (%ld) larger than M (%d)\n", (long) obs->get_num_symbols(), M)
5313  }
5314  }
5315 
5316  if (!reused_caches)
5317  {
5318 #ifdef USE_HMMPARALLEL_STRUCTURES
5319  for (int32_t i=0; i<parallel->get_num_threads(); i++)
5320  {
5321  SG_FREE(alpha_cache[i].table);
5322  SG_FREE(beta_cache[i].table);
5323  SG_FREE(states_per_observation_psi[i]);
5324  SG_FREE(path[i]);
5325 
5326  alpha_cache[i].table=NULL;
5327  beta_cache[i].table=NULL;
5328  states_per_observation_psi[i]=NULL;
5329  path[i]=NULL;
5330  } ;
5331 #else
5332  SG_FREE(alpha_cache.table);
5333  SG_FREE(beta_cache.table);
5334  SG_FREE(states_per_observation_psi);
5335  SG_FREE(path);
5336 
5337  alpha_cache.table=NULL;
5338  beta_cache.table=NULL;
5339  states_per_observation_psi=NULL;
5340  path=NULL;
5341 
5342 #endif //USE_HMMPARALLEL_STRUCTURES
5343  }
5344 
5345  if (obs!=NULL)
5346  {
5347  int32_t max_T=obs->get_max_vector_length();
5348 
5349  if (lambda)
5350  {
5351 #ifdef USE_HMMPARALLEL_STRUCTURES
5352  for (int32_t i=0; i<parallel->get_num_threads(); i++)
5353  {
5354  this->alpha_cache[i].table= lambda->alpha_cache[i].table;
5355  this->beta_cache[i].table= lambda->beta_cache[i].table;
5356  this->states_per_observation_psi[i]=lambda->states_per_observation_psi[i] ;
5357  this->path[i]=lambda->path[i];
5358  } ;
5359 #else
5360  this->alpha_cache.table= lambda->alpha_cache.table;
5361  this->beta_cache.table= lambda->beta_cache.table;
5362  this->states_per_observation_psi= lambda->states_per_observation_psi;
5363  this->path=lambda->path;
5364 #endif //USE_HMMPARALLEL_STRUCTURES
5365 
5366  this->reused_caches=true;
5367  }
5368  else
5369  {
5370  this->reused_caches=false;
5371 #ifdef USE_HMMPARALLEL_STRUCTURES
5372  SG_INFO("allocating mem for path-table of size %.2f Megabytes (%d*%d) each:\n", ((float32_t)max_T)*N*sizeof(T_STATES)/(1024*1024), max_T, N)
5373  for (int32_t i=0; i<parallel->get_num_threads(); i++)
5374  {
5375  if ((states_per_observation_psi[i]=SG_MALLOC(T_STATES,max_T*N))!=NULL)
5376  SG_DEBUG("path_table[%i] successfully allocated\n",i)
5377  else
5378  SG_ERROR("failed allocating memory for path_table[%i].\n",i)
5379  path[i]=SG_MALLOC(T_STATES, max_T);
5380  }
5381 #else // no USE_HMMPARALLEL_STRUCTURES
5382  SG_INFO("allocating mem of size %.2f Megabytes (%d*%d) for path-table ....", ((float32_t)max_T)*N*sizeof(T_STATES)/(1024*1024), max_T, N)
5383  if ((states_per_observation_psi=SG_MALLOC(T_STATES,max_T*N)) != NULL)
5384  SG_DONE()
5385  else
5386  SG_ERROR("failed.\n")
5387 
5388  path=SG_MALLOC(T_STATES, max_T);
5389 #endif // USE_HMMPARALLEL_STRUCTURES
5390 #ifdef USE_HMMCACHE
5391  SG_INFO("allocating mem for caches each of size %.2f Megabytes (%d*%d) ....\n", ((float32_t)max_T)*N*sizeof(T_ALPHA_BETA_TABLE)/(1024*1024), max_T, N)
5392 
5393 #ifdef USE_HMMPARALLEL_STRUCTURES
5394  for (int32_t i=0; i<parallel->get_num_threads(); i++)
5395  {
5396  if ((alpha_cache[i].table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N))!=NULL)
5397  SG_DEBUG("alpha_cache[%i].table successfully allocated\n",i)
5398  else
5399  SG_ERROR("allocation of alpha_cache[%i].table failed\n",i)
5400 
5401  if ((beta_cache[i].table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N)) != NULL)
5402  SG_DEBUG("beta_cache[%i].table successfully allocated\n",i)
5403  else
5404  SG_ERROR("allocation of beta_cache[%i].table failed\n",i)
5405  } ;
5406 #else // USE_HMMPARALLEL_STRUCTURES
5407  if ((alpha_cache.table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N)) != NULL)
5408  SG_DEBUG("alpha_cache.table successfully allocated\n")
5409  else
5410  SG_ERROR("allocation of alpha_cache.table failed\n")
5411 
5412  if ((beta_cache.table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N)) != NULL)
5413  SG_DEBUG("beta_cache.table successfully allocated\n")
5414  else
5415  SG_ERROR("allocation of beta_cache.table failed\n")
5416 
5417 #endif // USE_HMMPARALLEL_STRUCTURES
5418 #else // USE_HMMCACHE
5419 #ifdef USE_HMMPARALLEL_STRUCTURES
5420  for (int32_t i=0; i<parallel->get_num_threads(); i++)
5421  {
5422  alpha_cache[i].table=NULL ;
5423  beta_cache[i].table=NULL ;
5424  } ;
5425 #else //USE_HMMPARALLEL_STRUCTURES
5426  alpha_cache.table=NULL ;
5427  beta_cache.table=NULL ;
5428 #endif //USE_HMMPARALLEL_STRUCTURES
5429 #endif //USE_HMMCACHE
5430  }
5431  }
5432 
5433  //initialize pat/mod_prob as not calculated
5434  invalidate_model();
5435 }
5436 
5437 bool CHMM::permutation_entropy(int32_t window_width, int32_t sequence_number)
5438 {
5439  if (p_observations && window_width>0 &&
5440  ( sequence_number<0 || sequence_number < p_observations->get_num_vectors()))
5441  {
5442  int32_t min_sequence=sequence_number;
5443  int32_t max_sequence=sequence_number;
5444 
5445  if (sequence_number<0)
5446  {
5447  min_sequence=0;
5448  max_sequence=p_observations->get_num_vectors();
5449  SG_INFO("numseq: %d\n", max_sequence)
5450  }
5451 
5452  SG_INFO("min_sequence: %d max_sequence: %d\n", min_sequence, max_sequence)
5453  for (sequence_number=min_sequence; sequence_number<max_sequence; sequence_number++)
5454  {
5455  int32_t sequence_length=0;
5456  bool free_vec;
5457  uint16_t* obs=p_observations->get_feature_vector(sequence_number, sequence_length, free_vec);
5458 
5459  int32_t histsize=get_M();
5460  int64_t* hist=SG_MALLOC(int64_t, histsize);
5461  int32_t i,j;
5462 
5463  for (i=0; i<sequence_length-window_width; i++)
5464  {
5465  for (j=0; j<histsize; j++)
5466  hist[j]=0;
5467 
5468  uint16_t* ptr=&obs[i];
5469  for (j=0; j<window_width; j++)
5470  {
5471  hist[*ptr++]++;
5472  }
5473 
5474  float64_t perm_entropy=0;
5475  for (j=0; j<get_M(); j++)
5476  {
5477  float64_t p=
5478  (((float64_t) hist[j])+PSEUDO)/
5479  (window_width+get_M()*PSEUDO);
5480  perm_entropy+=p*log(p);
5481  }
5482 
5483  SG_PRINT("%f\n", perm_entropy)
5484  }
5485  p_observations->free_feature_vector(obs, sequence_number, free_vec);
5486 
5487  SG_FREE(hist);
5488  }
5489  return true;
5490  }
5491  else
5492  return false;
5493 }
5494 
5495 float64_t CHMM::get_log_derivative(int32_t num_param, int32_t num_example)
5496 {
5497  if (num_param<N)
5498  return model_derivative_p(num_param, num_example);
5499  else if (num_param<2*N)
5500  return model_derivative_q(num_param-N, num_example);
5501  else if (num_param<N*(N+2))
5502  {
5503  int32_t k=num_param-2*N;
5504  int32_t i=k/N;
5505  int32_t j=k%N;
5506  return model_derivative_a(i,j, num_example);
5507  }
5508  else if (num_param<N*(N+2+M))
5509  {
5510  int32_t k=num_param-N*(N+2);
5511  int32_t i=k/M;
5512  int32_t j=k%M;
5513  return model_derivative_b(i,j, num_example);
5514  }
5515 
5516  ASSERT(false)
5517  return -1;
5518 }
5519 
5520 float64_t CHMM::get_log_model_parameter(int32_t num_param)
5521 {
5522  if (num_param<N)
5523  return get_p(num_param);
5524  else if (num_param<2*N)
5525  return get_q(num_param-N);
5526  else if (num_param<N*(N+2))
5527  return transition_matrix_a[num_param-2*N];
5528  else if (num_param<N*(N+2+M))
5529  return observation_matrix_b[num_param-N*(N+2)];
5530 
5531  ASSERT(false)
5532  return -1;
5533 }
5534 
5535 
5536 //convergence criteria -tobeadjusted-
5537 bool CHMM::converged(float64_t x, float64_t y)
5538 {
5539  float64_t diff=y-x;
5540  float64_t absdiff=fabs(diff);
5541 
5542  SG_INFO("\n #%03d\tbest result so far: %G (eps: %f)", iteration_count, y, diff)
5543 
5544  if (iteration_count--==0 || (absdiff<epsilon && conv_it<=0))
5545  {
5546  iteration_count=iterations;
5547  SG_INFO("...finished\n")
5548  conv_it=5;
5549  return true;
5550  }
5551  else
5552  {
5553  if (absdiff<epsilon)
5554  conv_it--;
5555  else
5556  conv_it=5;
5557 
5558  return false;
5559  }
5560 }
5561 
5562 bool CHMM::baum_welch_viterbi_train(BaumWelchViterbiType type)
5563 {
5564  CHMM* estimate=new CHMM(this);
5565  CHMM* working=this;
5566  float64_t prob_max=-CMath::INFTY;
5567  float64_t prob=-CMath::INFTY;
5568  float64_t prob_train=CMath::ALMOST_NEG_INFTY;
5569  iteration_count=iterations;
5570 
5571  while (!converged(prob, prob_train) && (!CSignal::cancel_computations()))
5572  {
5573  CMath::swap(working, estimate);
5574  prob=prob_train;
5575 
5576  switch (type) {
5577  case BW_NORMAL:
5578  working->estimate_model_baum_welch(estimate); break;
5579  case BW_TRANS:
5580  working->estimate_model_baum_welch_trans(estimate); break;
5581  case BW_DEFINED:
5582  working->estimate_model_baum_welch_defined(estimate); break;
5583  case VIT_NORMAL:
5584  working->estimate_model_viterbi(estimate); break;
5585  case VIT_DEFINED:
5586  working->estimate_model_viterbi_defined(estimate); break;
5587  }
5588  prob_train=estimate->model_probability();
5589 
5590  if (prob_max<prob_train)
5591  prob_max=prob_train;
5592  }
5593 
5594  if (estimate == this)
5595  {
5596  estimate->copy_model(working);
5597  delete working;
5598  }
5599  else
5600  delete estimate;
5601 
5602  return true;
5603 }
SHOGUN Machine Learning Toolbox - Documentation