00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <shogun/distributions/HMM.h>
00012 #include <shogun/mathematics/Math.h>
00013 #include <shogun/io/SGIO.h>
00014 #include <shogun/lib/config.h>
00015 #include <shogun/lib/Signal.h>
00016 #include <shogun/base/Parallel.h>
00017 #include <shogun/features/StringFeatures.h>
00018 #include <shogun/features/Alphabet.h>
00019
00020 #include <stdlib.h>
00021 #include <stdio.h>
00022 #include <time.h>
00023 #include <ctype.h>
00024
00025 #define VAL_MACRO log((default_value == 0) ? (CMath::random(MIN_RAND, MAX_RAND)) : default_value)
00026 #define ARRAY_SIZE 65336
00027
00028 using namespace shogun;
00029
00031
00033
00034 const int32_t CHMM::GOTN= (1<<1);
00035 const int32_t CHMM::GOTM= (1<<2);
00036 const int32_t CHMM::GOTO= (1<<3);
00037 const int32_t CHMM::GOTa= (1<<4);
00038 const int32_t CHMM::GOTb= (1<<5);
00039 const int32_t CHMM::GOTp= (1<<6);
00040 const int32_t CHMM::GOTq= (1<<7);
00041
00042 const int32_t CHMM::GOTlearn_a= (1<<1);
00043 const int32_t CHMM::GOTlearn_b= (1<<2);
00044 const int32_t CHMM::GOTlearn_p= (1<<3);
00045 const int32_t CHMM::GOTlearn_q= (1<<4);
00046 const int32_t CHMM::GOTconst_a= (1<<5);
00047 const int32_t CHMM::GOTconst_b= (1<<6);
00048 const int32_t CHMM::GOTconst_p= (1<<7);
00049 const int32_t CHMM::GOTconst_q= (1<<8);
00050
00051 enum E_STATE
00052 {
00053 INITIAL,
00054 ARRAYs,
00055 GET_N,
00056 GET_M,
00057 GET_a,
00058 GET_b,
00059 GET_p,
00060 GET_q,
00061 GET_learn_a,
00062 GET_learn_b,
00063 GET_learn_p,
00064 GET_learn_q,
00065 GET_const_a,
00066 GET_const_b,
00067 GET_const_p,
00068 GET_const_q,
00069 COMMENT,
00070 END
00071 };
00072
00073
00074 #ifdef FIX_POS
00075 const char Model::FIX_DISALLOWED=0 ;
00076 const char Model::FIX_ALLOWED=1 ;
00077 const char Model::FIX_DEFAULT=-1 ;
00078 const float64_t Model::DISALLOWED_PENALTY=CMath::ALMOST_NEG_INFTY ;
00079 #endif
00080
00081 Model::Model()
00082 {
00083 const_a=SG_MALLOC(int, ARRAY_SIZE);
00084 const_b=SG_MALLOC(int, ARRAY_SIZE);
00085 const_p=SG_MALLOC(int, ARRAY_SIZE);
00086 const_q=SG_MALLOC(int, ARRAY_SIZE);
00087 const_a_val=SG_MALLOC(float64_t, ARRAY_SIZE);
00088 const_b_val=SG_MALLOC(float64_t, ARRAY_SIZE);
00089 const_p_val=SG_MALLOC(float64_t, ARRAY_SIZE);
00090 const_q_val=SG_MALLOC(float64_t, ARRAY_SIZE);
00091
00092
00093 learn_a=SG_MALLOC(int, ARRAY_SIZE);
00094 learn_b=SG_MALLOC(int, ARRAY_SIZE);
00095 learn_p=SG_MALLOC(int, ARRAY_SIZE);
00096 learn_q=SG_MALLOC(int, ARRAY_SIZE);
00097
00098 #ifdef FIX_POS
00099 fix_pos_state = SG_MALLOC(char, ARRAY_SIZE);
00100 #endif
00101 for (int32_t i=0; i<ARRAY_SIZE; i++)
00102 {
00103 const_a[i]=-1 ;
00104 const_b[i]=-1 ;
00105 const_p[i]=-1 ;
00106 const_q[i]=-1 ;
00107 const_a_val[i]=1.0 ;
00108 const_b_val[i]=1.0 ;
00109 const_p_val[i]=1.0 ;
00110 const_q_val[i]=1.0 ;
00111 learn_a[i]=-1 ;
00112 learn_b[i]=-1 ;
00113 learn_p[i]=-1 ;
00114 learn_q[i]=-1 ;
00115 #ifdef FIX_POS
00116 fix_pos_state[i] = FIX_DEFAULT ;
00117 #endif
00118 } ;
00119 }
00120
00121 Model::~Model()
00122 {
00123 SG_FREE(const_a);
00124 SG_FREE(const_b);
00125 SG_FREE(const_p);
00126 SG_FREE(const_q);
00127 SG_FREE(const_a_val);
00128 SG_FREE(const_b_val);
00129 SG_FREE(const_p_val);
00130 SG_FREE(const_q_val);
00131
00132 SG_FREE(learn_a);
00133 SG_FREE(learn_b);
00134 SG_FREE(learn_p);
00135 SG_FREE(learn_q);
00136
00137 #ifdef FIX_POS
00138 SG_FREE(fix_pos_state);
00139 #endif
00140
00141 }
00142
00143 CHMM::CHMM()
00144 {
00145 N=0;
00146 M=0;
00147 model=NULL;
00148 status=false;
00149 }
00150
00151 CHMM::CHMM(CHMM* h)
00152 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
00153 {
00154 #ifdef USE_HMMPARALLEL_STRUCTURES
00155 SG_INFO( "hmm is using %i separate tables\n", parallel->get_num_threads()) ;
00156 #endif
00157
00158 this->N=h->get_N();
00159 this->M=h->get_M();
00160 status=initialize(NULL, h->get_pseudo());
00161 this->copy_model(h);
00162 set_observations(h->p_observations);
00163 }
00164
00165 CHMM::CHMM(int32_t p_N, int32_t p_M, Model* p_model, float64_t p_PSEUDO)
00166 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
00167 {
00168 this->N=p_N;
00169 this->M=p_M;
00170 model=NULL ;
00171
00172 #ifdef USE_HMMPARALLEL_STRUCTURES
00173 SG_INFO( "hmm is using %i separate tables\n", parallel->get_num_threads()) ;
00174 #endif
00175
00176 status=initialize(p_model, p_PSEUDO);
00177 }
00178
00179 CHMM::CHMM(
00180 CStringFeatures<uint16_t>* obs, int32_t p_N, int32_t p_M,
00181 float64_t p_PSEUDO)
00182 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
00183 {
00184 this->N=p_N;
00185 this->M=p_M;
00186 model=NULL ;
00187
00188 #ifdef USE_HMMPARALLEL_STRUCTURES
00189 SG_INFO( "hmm is using %i separate tables\n", parallel->get_num_threads()) ;
00190 #endif
00191
00192 initialize(model, p_PSEUDO);
00193 set_observations(obs);
00194 }
00195
00196 CHMM::CHMM(int32_t p_N, float64_t* p, float64_t* q, float64_t* a)
00197 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
00198 {
00199 this->N=p_N;
00200 this->M=0;
00201 model=NULL ;
00202
00203 trans_list_forward = NULL ;
00204 trans_list_forward_cnt = NULL ;
00205 trans_list_forward_val = NULL ;
00206 trans_list_backward = NULL ;
00207 trans_list_backward_cnt = NULL ;
00208 trans_list_len = 0 ;
00209 mem_initialized = false ;
00210
00211 this->transition_matrix_a=NULL;
00212 this->observation_matrix_b=NULL;
00213 this->initial_state_distribution_p=NULL;
00214 this->end_state_distribution_q=NULL;
00215 this->PSEUDO= PSEUDO;
00216 this->model= model;
00217 this->p_observations=NULL;
00218 this->reused_caches=false;
00219
00220 #ifdef USE_HMMPARALLEL_STRUCTURES
00221 this->alpha_cache=NULL;
00222 this->beta_cache=NULL;
00223 #else
00224 this->alpha_cache.table=NULL;
00225 this->beta_cache.table=NULL;
00226 this->alpha_cache.dimension=0;
00227 this->beta_cache.dimension=0;
00228 #endif
00229
00230 this->states_per_observation_psi=NULL ;
00231 this->path=NULL;
00232 arrayN1=NULL ;
00233 arrayN2=NULL ;
00234
00235 this->loglikelihood=false;
00236 mem_initialized = true ;
00237
00238 transition_matrix_a=a ;
00239 observation_matrix_b=NULL ;
00240 initial_state_distribution_p=p ;
00241 end_state_distribution_q=q ;
00242 transition_matrix_A=NULL ;
00243 observation_matrix_B=NULL ;
00244
00245
00246 }
00247
00248 CHMM::CHMM(
00249 int32_t p_N, float64_t* p, float64_t* q, int32_t num_trans,
00250 float64_t* a_trans)
00251 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
00252 {
00253 model=NULL ;
00254
00255 this->N=p_N;
00256 this->M=0;
00257
00258 trans_list_forward = NULL ;
00259 trans_list_forward_cnt = NULL ;
00260 trans_list_forward_val = NULL ;
00261 trans_list_backward = NULL ;
00262 trans_list_backward_cnt = NULL ;
00263 trans_list_len = 0 ;
00264 mem_initialized = false ;
00265
00266 this->transition_matrix_a=NULL;
00267 this->observation_matrix_b=NULL;
00268 this->initial_state_distribution_p=NULL;
00269 this->end_state_distribution_q=NULL;
00270 this->PSEUDO= PSEUDO;
00271 this->model= model;
00272 this->p_observations=NULL;
00273 this->reused_caches=false;
00274
00275 #ifdef USE_HMMPARALLEL_STRUCTURES
00276 this->alpha_cache=NULL;
00277 this->beta_cache=NULL;
00278 #else
00279 this->alpha_cache.table=NULL;
00280 this->beta_cache.table=NULL;
00281 this->alpha_cache.dimension=0;
00282 this->beta_cache.dimension=0;
00283 #endif
00284
00285 this->states_per_observation_psi=NULL ;
00286 this->path=NULL;
00287 arrayN1=NULL ;
00288 arrayN2=NULL ;
00289
00290 this->loglikelihood=false;
00291 mem_initialized = true ;
00292
00293 trans_list_forward_cnt=NULL ;
00294 trans_list_len = N ;
00295 trans_list_forward = SG_MALLOC(T_STATES*, N);
00296 trans_list_forward_val = SG_MALLOC(float64_t*, N);
00297 trans_list_forward_cnt = SG_MALLOC(T_STATES, N);
00298
00299 int32_t start_idx=0;
00300 for (int32_t j=0; j<N; j++)
00301 {
00302 int32_t old_start_idx=start_idx;
00303
00304 while (start_idx<num_trans && a_trans[start_idx+num_trans]==j)
00305 {
00306 start_idx++;
00307
00308 if (start_idx>1 && start_idx<num_trans)
00309 ASSERT(a_trans[start_idx+num_trans-1]<=
00310 a_trans[start_idx+num_trans]);
00311 }
00312
00313 if (start_idx>1 && start_idx<num_trans)
00314 ASSERT(a_trans[start_idx+num_trans-1]<=
00315 a_trans[start_idx+num_trans]);
00316
00317 int32_t len=start_idx-old_start_idx;
00318 ASSERT(len>=0);
00319
00320 trans_list_forward_cnt[j] = 0 ;
00321
00322 if (len>0)
00323 {
00324 trans_list_forward[j] = SG_MALLOC(T_STATES, len);
00325 trans_list_forward_val[j] = SG_MALLOC(float64_t, len);
00326 }
00327 else
00328 {
00329 trans_list_forward[j] = NULL;
00330 trans_list_forward_val[j] = NULL;
00331 }
00332 }
00333
00334 for (int32_t i=0; i<num_trans; i++)
00335 {
00336 int32_t from = (int32_t)a_trans[i+num_trans] ;
00337 int32_t to = (int32_t)a_trans[i] ;
00338 float64_t val = a_trans[i+num_trans*2] ;
00339
00340 ASSERT(from>=0 && from<N);
00341 ASSERT(to>=0 && to<N);
00342
00343 trans_list_forward[from][trans_list_forward_cnt[from]]=to ;
00344 trans_list_forward_val[from][trans_list_forward_cnt[from]]=val ;
00345 trans_list_forward_cnt[from]++ ;
00346
00347 } ;
00348
00349 transition_matrix_a=NULL ;
00350 observation_matrix_b=NULL ;
00351 initial_state_distribution_p=p ;
00352 end_state_distribution_q=q ;
00353 transition_matrix_A=NULL ;
00354 observation_matrix_B=NULL ;
00355
00356
00357 }
00358
00359
00360 CHMM::CHMM(FILE* model_file, float64_t p_PSEUDO)
00361 : CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
00362 {
00363 #ifdef USE_HMMPARALLEL_STRUCTURES
00364 SG_INFO( "hmm is using %i separate tables\n", parallel->get_num_threads()) ;
00365 #endif
00366
00367 status=initialize(NULL, p_PSEUDO, model_file);
00368 }
00369
00370 CHMM::~CHMM()
00371 {
00372 SG_UNREF(p_observations);
00373
00374 if (trans_list_forward_cnt)
00375 SG_FREE(trans_list_forward_cnt);
00376 if (trans_list_backward_cnt)
00377 SG_FREE(trans_list_backward_cnt);
00378 if (trans_list_forward)
00379 {
00380 for (int32_t i=0; i<trans_list_len; i++)
00381 if (trans_list_forward[i])
00382 SG_FREE(trans_list_forward[i]);
00383 SG_FREE(trans_list_forward);
00384 }
00385 if (trans_list_forward_val)
00386 {
00387 for (int32_t i=0; i<trans_list_len; i++)
00388 if (trans_list_forward_val[i])
00389 SG_FREE(trans_list_forward_val[i]);
00390 SG_FREE(trans_list_forward_val);
00391 }
00392 if (trans_list_backward)
00393 {
00394 for (int32_t i=0; i<trans_list_len; i++)
00395 if (trans_list_backward[i])
00396 SG_FREE(trans_list_backward[i]);
00397 SG_FREE(trans_list_backward);
00398 } ;
00399
00400 free_state_dependend_arrays();
00401
00402 if (!reused_caches)
00403 {
00404 #ifdef USE_HMMPARALLEL_STRUCTURES
00405 for (int32_t i=0; i<parallel->get_num_threads(); i++)
00406 {
00407 SG_FREE(alpha_cache[i].table);
00408 SG_FREE(beta_cache[i].table);
00409 alpha_cache[i].table=NULL;
00410 beta_cache[i].table=NULL;
00411 }
00412
00413 SG_FREE(alpha_cache);
00414 SG_FREE(beta_cache);
00415 alpha_cache=NULL;
00416 beta_cache=NULL;
00417 #else // USE_HMMPARALLEL_STRUCTURES
00418 SG_FREE(alpha_cache.table);
00419 SG_FREE(beta_cache.table);
00420 alpha_cache.table=NULL;
00421 beta_cache.table=NULL;
00422 #endif // USE_HMMPARALLEL_STRUCTURES
00423
00424 SG_FREE(states_per_observation_psi);
00425 states_per_observation_psi=NULL;
00426 }
00427
00428 #ifdef USE_LOGSUMARRAY
00429 #ifdef USE_HMMPARALLEL_STRUCTURES
00430 {
00431 for (int32_t i=0; i<parallel->get_num_threads(); i++)
00432 SG_FREE(arrayS[i]);
00433 SG_FREE(arrayS);
00434 } ;
00435 #else //USE_HMMPARALLEL_STRUCTURES
00436 SG_FREE(arrayS);
00437 #endif //USE_HMMPARALLEL_STRUCTURES
00438 #endif //USE_LOGSUMARRAY
00439
00440 if (!reused_caches)
00441 {
00442 #ifdef USE_HMMPARALLEL_STRUCTURES
00443 SG_FREE(path_prob_updated);
00444 SG_FREE(path_prob_dimension);
00445 for (int32_t i=0; i<parallel->get_num_threads(); i++)
00446 SG_FREE(path[i]);
00447 #endif //USE_HMMPARALLEL_STRUCTURES
00448 SG_FREE(path);
00449 }
00450 }
00451
00452 bool CHMM::train(CFeatures* data)
00453 {
00454 if (data)
00455 {
00456 if (data->get_feature_class() != C_STRING ||
00457 data->get_feature_type() != F_WORD)
00458 {
00459 SG_ERROR("Expected features of class string type word\n");
00460 }
00461 set_observations((CStringFeatures<uint16_t>*) data);
00462 }
00463 return baum_welch_viterbi_train(BW_NORMAL);
00464 }
00465
00466 bool CHMM::alloc_state_dependend_arrays()
00467 {
00468
00469 if (!transition_matrix_a && !observation_matrix_b &&
00470 !initial_state_distribution_p && !end_state_distribution_q)
00471 {
00472 transition_matrix_a=SG_MALLOC(float64_t, N*N);
00473 observation_matrix_b=SG_MALLOC(float64_t, N*M);
00474 initial_state_distribution_p=SG_MALLOC(float64_t, N);
00475 end_state_distribution_q=SG_MALLOC(float64_t, N);
00476 init_model_random();
00477 convert_to_log();
00478 }
00479
00480 #ifdef USE_HMMPARALLEL_STRUCTURES
00481 for (int32_t i=0; i<parallel->get_num_threads(); i++)
00482 {
00483 arrayN1[i]=SG_MALLOC(float64_t, N);
00484 arrayN2[i]=SG_MALLOC(float64_t, N);
00485 }
00486 #else //USE_HMMPARALLEL_STRUCTURES
00487 arrayN1=SG_MALLOC(float64_t, N);
00488 arrayN2=SG_MALLOC(float64_t, N);
00489 #endif //USE_HMMPARALLEL_STRUCTURES
00490
00491 #ifdef LOG_SUMARRAY
00492 #ifdef USE_HMMPARALLEL_STRUCTURES
00493 for (int32_t i=0; i<parallel->get_num_threads(); i++)
00494 arrayS[i]=SG_MALLOC(float64_t, (int32_t)(this->N/2+1));
00495 #else //USE_HMMPARALLEL_STRUCTURES
00496 arrayS=SG_MALLOC(float64_t, (int32_t)(this->N/2+1));
00497 #endif //USE_HMMPARALLEL_STRUCTURES
00498 #endif //LOG_SUMARRAY
00499 transition_matrix_A=SG_MALLOC(float64_t, this->N*this->N);
00500 observation_matrix_B=SG_MALLOC(float64_t, this->N*this->M);
00501
00502 if (p_observations)
00503 {
00504 #ifdef USE_HMMPARALLEL_STRUCTURES
00505 if (alpha_cache[0].table!=NULL)
00506 #else //USE_HMMPARALLEL_STRUCTURES
00507 if (alpha_cache.table!=NULL)
00508 #endif //USE_HMMPARALLEL_STRUCTURES
00509 set_observations(p_observations);
00510 else
00511 set_observation_nocache(p_observations);
00512 SG_UNREF(p_observations);
00513 }
00514
00515 this->invalidate_model();
00516
00517 return ((transition_matrix_A != NULL) && (observation_matrix_B != NULL) &&
00518 (transition_matrix_a != NULL) && (observation_matrix_b != NULL) &&
00519 (initial_state_distribution_p != NULL) &&
00520 (end_state_distribution_q != NULL));
00521 }
00522
00523 void CHMM::free_state_dependend_arrays()
00524 {
00525 #ifdef USE_HMMPARALLEL_STRUCTURES
00526 for (int32_t i=0; i<parallel->get_num_threads(); i++)
00527 {
00528 SG_FREE(arrayN1[i]);
00529 SG_FREE(arrayN2[i]);
00530
00531 arrayN1[i]=NULL;
00532 arrayN2[i]=NULL;
00533 }
00534 #else
00535 SG_FREE(arrayN1);
00536 SG_FREE(arrayN2);
00537 arrayN1=NULL;
00538 arrayN2=NULL;
00539 #endif
00540 if (observation_matrix_b)
00541 {
00542 SG_FREE(transition_matrix_A);
00543 SG_FREE(observation_matrix_B);
00544 SG_FREE(transition_matrix_a);
00545 SG_FREE(observation_matrix_b);
00546 SG_FREE(initial_state_distribution_p);
00547 SG_FREE(end_state_distribution_q);
00548 } ;
00549
00550 transition_matrix_A=NULL;
00551 observation_matrix_B=NULL;
00552 transition_matrix_a=NULL;
00553 observation_matrix_b=NULL;
00554 initial_state_distribution_p=NULL;
00555 end_state_distribution_q=NULL;
00556 }
00557
00558 bool CHMM::initialize(Model* m, float64_t pseudo, FILE* modelfile)
00559 {
00560
00561 bool files_ok=true;
00562
00563 trans_list_forward = NULL ;
00564 trans_list_forward_cnt = NULL ;
00565 trans_list_forward_val = NULL ;
00566 trans_list_backward = NULL ;
00567 trans_list_backward_cnt = NULL ;
00568 trans_list_len = 0 ;
00569 mem_initialized = false ;
00570
00571 this->transition_matrix_a=NULL;
00572 this->observation_matrix_b=NULL;
00573 this->initial_state_distribution_p=NULL;
00574 this->end_state_distribution_q=NULL;
00575 this->PSEUDO= pseudo;
00576 this->model= m;
00577 this->p_observations=NULL;
00578 this->reused_caches=false;
00579
00580 #ifdef USE_HMMPARALLEL_STRUCTURES
00581 alpha_cache=SG_MALLOC(T_ALPHA_BETA, parallel->get_num_threads());
00582 beta_cache=SG_MALLOC(T_ALPHA_BETA, parallel->get_num_threads());
00583 states_per_observation_psi=SG_MALLOC(P_STATES, parallel->get_num_threads());
00584
00585 for (int32_t i=0; i<parallel->get_num_threads(); i++)
00586 {
00587 this->alpha_cache[i].table=NULL;
00588 this->beta_cache[i].table=NULL;
00589 this->alpha_cache[i].dimension=0;
00590 this->beta_cache[i].dimension=0;
00591 this->states_per_observation_psi[i]=NULL ;
00592 }
00593
00594 #else // USE_HMMPARALLEL_STRUCTURES
00595 this->alpha_cache.table=NULL;
00596 this->beta_cache.table=NULL;
00597 this->alpha_cache.dimension=0;
00598 this->beta_cache.dimension=0;
00599 this->states_per_observation_psi=NULL ;
00600 #endif //USE_HMMPARALLEL_STRUCTURES
00601
00602 if (modelfile)
00603 files_ok= files_ok && load_model(modelfile);
00604
00605 #ifdef USE_HMMPARALLEL_STRUCTURES
00606 path_prob_updated=SG_MALLOC(bool, parallel->get_num_threads());
00607 path_prob_dimension=SG_MALLOC(int, parallel->get_num_threads());
00608
00609 path=SG_MALLOC(P_STATES, parallel->get_num_threads());
00610
00611 for (int32_t i=0; i<parallel->get_num_threads(); i++)
00612 this->path[i]=NULL;
00613
00614 #else // USE_HMMPARALLEL_STRUCTURES
00615 this->path=NULL;
00616
00617 #endif //USE_HMMPARALLEL_STRUCTURES
00618
00619 #ifdef USE_HMMPARALLEL_STRUCTURES
00620 arrayN1=SG_MALLOC(float64_t*, parallel->get_num_threads());
00621 arrayN2=SG_MALLOC(float64_t*, parallel->get_num_threads());
00622 #endif //USE_HMMPARALLEL_STRUCTURES
00623
00624 #ifdef LOG_SUMARRAY
00625 #ifdef USE_HMMPARALLEL_STRUCTURES
00626 arrayS=SG_MALLOC(float64_t*, parallel->get_num_threads());
00627 #endif // USE_HMMPARALLEL_STRUCTURES
00628 #endif //LOG_SUMARRAY
00629
00630 alloc_state_dependend_arrays();
00631
00632 this->loglikelihood=false;
00633 mem_initialized = true ;
00634 this->invalidate_model();
00635
00636 return ((files_ok) &&
00637 (transition_matrix_A != NULL) && (observation_matrix_B != NULL) &&
00638 (transition_matrix_a != NULL) && (observation_matrix_b != NULL) && (initial_state_distribution_p != NULL) &&
00639 (end_state_distribution_q != NULL));
00640 }
00641
00642
00643
00644
00645
00646
00647 float64_t CHMM::forward_comp(int32_t time, int32_t state, int32_t dimension)
00648 {
00649 T_ALPHA_BETA_TABLE* alpha_new;
00650 T_ALPHA_BETA_TABLE* alpha;
00651 T_ALPHA_BETA_TABLE* dummy;
00652 if (time<1)
00653 time=0;
00654
00655
00656 int32_t wanted_time=time;
00657
00658 if (ALPHA_CACHE(dimension).table)
00659 {
00660 alpha=&ALPHA_CACHE(dimension).table[0];
00661 alpha_new=&ALPHA_CACHE(dimension).table[N];
00662 time=p_observations->get_vector_length(dimension)+1;
00663 }
00664 else
00665 {
00666 alpha_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
00667 alpha=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
00668 }
00669
00670 if (time<1)
00671 return get_p(state) + get_b(state, p_observations->get_feature(dimension,0));
00672 else
00673 {
00674
00675 for (int32_t i=0; i<N; i++)
00676 alpha[i] = get_p(i) + get_b(i, p_observations->get_feature(dimension,0)) ;
00677
00678
00679 for (register int32_t t=1; t<time && t < p_observations->get_vector_length(dimension); t++)
00680 {
00681
00682 for (int32_t j=0; j<N; j++)
00683 {
00684 register int32_t i, num = trans_list_forward_cnt[j] ;
00685 float64_t sum=-CMath::INFTY;
00686 for (i=0; i < num; i++)
00687 {
00688 int32_t ii = trans_list_forward[j][i] ;
00689 sum = CMath::logarithmic_sum(sum, alpha[ii] + get_a(ii,j));
00690 } ;
00691
00692 alpha_new[j]= sum + get_b(j, p_observations->get_feature(dimension,t));
00693 }
00694
00695 if (!ALPHA_CACHE(dimension).table)
00696 {
00697 dummy=alpha;
00698 alpha=alpha_new;
00699 alpha_new=dummy;
00700 }
00701 else
00702 {
00703 alpha=alpha_new;
00704 alpha_new+=N;
00705 }
00706 }
00707
00708
00709 if (time<p_observations->get_vector_length(dimension))
00710 {
00711 register int32_t i, num=trans_list_forward_cnt[state];
00712 register float64_t sum=-CMath::INFTY;
00713 for (i=0; i<num; i++)
00714 {
00715 int32_t ii = trans_list_forward[state][i] ;
00716 sum= CMath::logarithmic_sum(sum, alpha[ii] + get_a(ii, state));
00717 } ;
00718
00719 return sum + get_b(state, p_observations->get_feature(dimension,time));
00720 }
00721 else
00722 {
00723
00724 register int32_t i ;
00725 float64_t sum ;
00726 sum=-CMath::INFTY;
00727 for (i=0; i<N; i++)
00728 sum=CMath::logarithmic_sum(sum, alpha[i] + get_q(i));
00729
00730 if (!ALPHA_CACHE(dimension).table)
00731 return sum;
00732 else
00733 {
00734 ALPHA_CACHE(dimension).dimension=dimension;
00735 ALPHA_CACHE(dimension).updated=true;
00736 ALPHA_CACHE(dimension).sum=sum;
00737
00738 if (wanted_time<p_observations->get_vector_length(dimension))
00739 return ALPHA_CACHE(dimension).table[wanted_time*N+state];
00740 else
00741 return ALPHA_CACHE(dimension).sum;
00742 }
00743 }
00744 }
00745 }
00746
00747
00748
00749
00750
00751 float64_t CHMM::forward_comp_old(int32_t time, int32_t state, int32_t dimension)
00752 {
00753 T_ALPHA_BETA_TABLE* alpha_new;
00754 T_ALPHA_BETA_TABLE* alpha;
00755 T_ALPHA_BETA_TABLE* dummy;
00756 if (time<1)
00757 time=0;
00758
00759 int32_t wanted_time=time;
00760
00761 if (ALPHA_CACHE(dimension).table)
00762 {
00763 alpha=&ALPHA_CACHE(dimension).table[0];
00764 alpha_new=&ALPHA_CACHE(dimension).table[N];
00765 time=p_observations->get_vector_length(dimension)+1;
00766 }
00767 else
00768 {
00769 alpha_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
00770 alpha=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
00771 }
00772
00773 if (time<1)
00774 return get_p(state) + get_b(state, p_observations->get_feature(dimension,0));
00775 else
00776 {
00777
00778 for (int32_t i=0; i<N; i++)
00779 alpha[i] = get_p(i) + get_b(i, p_observations->get_feature(dimension,0)) ;
00780
00781
00782 for (register int32_t t=1; t<time && t < p_observations->get_vector_length(dimension); t++)
00783 {
00784
00785 for (int32_t j=0; j<N; j++)
00786 {
00787 register int32_t i ;
00788 #ifdef USE_LOGSUMARRAY
00789 for (i=0; i<(N>>1); i++)
00790 ARRAYS(dimension)[i]=CMath::logarithmic_sum(alpha[i<<1] + get_a(i<<1,j),
00791 alpha[(i<<1)+1] + get_a((i<<1)+1,j));
00792 if (N%2==1)
00793 alpha_new[j]=get_b(j, p_observations->get_feature(dimension,t))+
00794 CMath::logarithmic_sum(alpha[N-1]+get_a(N-1,j),
00795 CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
00796 else
00797 alpha_new[j]=get_b(j, p_observations->get_feature(dimension,t))+CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
00798 #else //USE_LOGSUMARRAY
00799 float64_t sum=-CMath::INFTY;
00800 for (i=0; i<N; i++)
00801 sum= CMath::logarithmic_sum(sum, alpha[i] + get_a(i,j));
00802
00803 alpha_new[j]= sum + get_b(j, p_observations->get_feature(dimension,t));
00804 #endif //USE_LOGSUMARRAY
00805 }
00806
00807 if (!ALPHA_CACHE(dimension).table)
00808 {
00809 dummy=alpha;
00810 alpha=alpha_new;
00811 alpha_new=dummy;
00812 }
00813 else
00814 {
00815 alpha=alpha_new;
00816 alpha_new+=N;
00817 }
00818 }
00819
00820
00821 if (time<p_observations->get_vector_length(dimension))
00822 {
00823 register int32_t i;
00824 #ifdef USE_LOGSUMARRAY
00825 for (i=0; i<(N>>1); i++)
00826 ARRAYS(dimension)[i]=CMath::logarithmic_sum(alpha[i<<1] + get_a(i<<1,state),
00827 alpha[(i<<1)+1] + get_a((i<<1)+1,state));
00828 if (N%2==1)
00829 return get_b(state, p_observations->get_feature(dimension,time))+
00830 CMath::logarithmic_sum(alpha[N-1]+get_a(N-1,state),
00831 CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
00832 else
00833 return get_b(state, p_observations->get_feature(dimension,time))+CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
00834 #else //USE_LOGSUMARRAY
00835 register float64_t sum=-CMath::INFTY;
00836 for (i=0; i<N; i++)
00837 sum= CMath::logarithmic_sum(sum, alpha[i] + get_a(i, state));
00838
00839 return sum + get_b(state, p_observations->get_feature(dimension,time));
00840 #endif //USE_LOGSUMARRAY
00841 }
00842 else
00843 {
00844
00845 register int32_t i ;
00846 float64_t sum ;
00847 #ifdef USE_LOGSUMARRAY
00848 for (i=0; i<(N>>1); i++)
00849 ARRAYS(dimension)[i]=CMath::logarithmic_sum(alpha[i<<1] + get_q(i<<1),
00850 alpha[(i<<1)+1] + get_q((i<<1)+1));
00851 if (N%2==1)
00852 sum=CMath::logarithmic_sum(alpha[N-1]+get_q(N-1),
00853 CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
00854 else
00855 sum=CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
00856 #else //USE_LOGSUMARRAY
00857 sum=-CMath::INFTY;
00858 for (i=0; i<N; i++)
00859 sum=CMath::logarithmic_sum(sum, alpha[i] + get_q(i));
00860 #endif //USE_LOGSUMARRAY
00861
00862 if (!ALPHA_CACHE(dimension).table)
00863 return sum;
00864 else
00865 {
00866 ALPHA_CACHE(dimension).dimension=dimension;
00867 ALPHA_CACHE(dimension).updated=true;
00868 ALPHA_CACHE(dimension).sum=sum;
00869
00870 if (wanted_time<p_observations->get_vector_length(dimension))
00871 return ALPHA_CACHE(dimension).table[wanted_time*N+state];
00872 else
00873 return ALPHA_CACHE(dimension).sum;
00874 }
00875 }
00876 }
00877 }
00878
00879
00880
00881
00882
00883 float64_t CHMM::backward_comp(int32_t time, int32_t state, int32_t dimension)
00884 {
00885 T_ALPHA_BETA_TABLE* beta_new;
00886 T_ALPHA_BETA_TABLE* beta;
00887 T_ALPHA_BETA_TABLE* dummy;
00888 int32_t wanted_time=time;
00889
00890 if (time<0)
00891 forward(time, state, dimension);
00892
00893 if (BETA_CACHE(dimension).table)
00894 {
00895 beta=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-1)];
00896 beta_new=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-2)];
00897 time=-1;
00898 }
00899 else
00900 {
00901 beta_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
00902 beta=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
00903 }
00904
00905 if (time>=p_observations->get_vector_length(dimension)-1)
00906
00907
00908 return get_q(state);
00909 else
00910 {
00911
00912 for (register int32_t i=0; i<N; i++)
00913 beta[i]=get_q(i);
00914
00915
00916 for (register int32_t t=p_observations->get_vector_length(dimension)-1; t>time+1 && t>0; t--)
00917 {
00918 for (register int32_t i=0; i<N; i++)
00919 {
00920 register int32_t j, num=trans_list_backward_cnt[i] ;
00921 float64_t sum=-CMath::INFTY;
00922 for (j=0; j<num; j++)
00923 {
00924 int32_t jj = trans_list_backward[i][j] ;
00925 sum= CMath::logarithmic_sum(sum, get_a(i, jj) + get_b(jj, p_observations->get_feature(dimension,t)) + beta[jj]);
00926 } ;
00927 beta_new[i]=sum;
00928 }
00929
00930 if (!BETA_CACHE(dimension).table)
00931 {
00932 dummy=beta;
00933 beta=beta_new;
00934 beta_new=dummy;
00935 }
00936 else
00937 {
00938 beta=beta_new;
00939 beta_new-=N;
00940 }
00941 }
00942
00943 if (time>=0)
00944 {
00945 register int32_t j, num=trans_list_backward_cnt[state] ;
00946 float64_t sum=-CMath::INFTY;
00947 for (j=0; j<num; j++)
00948 {
00949 int32_t jj = trans_list_backward[state][j] ;
00950 sum= CMath::logarithmic_sum(sum, get_a(state, jj) + get_b(jj, p_observations->get_feature(dimension,time+1))+beta[jj]);
00951 } ;
00952 return sum;
00953 }
00954 else
00955 {
00956 if (BETA_CACHE(dimension).table)
00957 {
00958 float64_t sum=-CMath::INFTY;
00959 for (register int32_t j=0; j<N; j++)
00960 sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
00961 BETA_CACHE(dimension).sum=sum;
00962 BETA_CACHE(dimension).dimension=dimension;
00963 BETA_CACHE(dimension).updated=true;
00964
00965 if (wanted_time<p_observations->get_vector_length(dimension))
00966 return BETA_CACHE(dimension).table[wanted_time*N+state];
00967 else
00968 return BETA_CACHE(dimension).sum;
00969 }
00970 else
00971 {
00972 float64_t sum=-CMath::INFTY;
00973 for (register int32_t j=0; j<N; j++)
00974 sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
00975 return sum;
00976 }
00977 }
00978 }
00979 }
00980
00981
00982 float64_t CHMM::backward_comp_old(
00983 int32_t time, int32_t state, int32_t dimension)
00984 {
00985 T_ALPHA_BETA_TABLE* beta_new;
00986 T_ALPHA_BETA_TABLE* beta;
00987 T_ALPHA_BETA_TABLE* dummy;
00988 int32_t wanted_time=time;
00989
00990 if (time<0)
00991 forward(time, state, dimension);
00992
00993 if (BETA_CACHE(dimension).table)
00994 {
00995 beta=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-1)];
00996 beta_new=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-2)];
00997 time=-1;
00998 }
00999 else
01000 {
01001 beta_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
01002 beta=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
01003 }
01004
01005 if (time>=p_observations->get_vector_length(dimension)-1)
01006
01007
01008 return get_q(state);
01009 else
01010 {
01011
01012 for (register int32_t i=0; i<N; i++)
01013 beta[i]=get_q(i);
01014
01015
01016 for (register int32_t t=p_observations->get_vector_length(dimension)-1; t>time+1 && t>0; t--)
01017 {
01018 for (register int32_t i=0; i<N; i++)
01019 {
01020 register int32_t j ;
01021 #ifdef USE_LOGSUMARRAY
01022 for (j=0; j<(N>>1); j++)
01023 ARRAYS(dimension)[j]=CMath::logarithmic_sum(
01024 get_a(i, j<<1) + get_b(j<<1, p_observations->get_feature(dimension,t)) + beta[j<<1],
01025 get_a(i, (j<<1)+1) + get_b((j<<1)+1, p_observations->get_feature(dimension,t)) + beta[(j<<1)+1]);
01026 if (N%2==1)
01027 beta_new[i]=CMath::logarithmic_sum(get_a(i, N-1) + get_b(N-1, p_observations->get_feature(dimension,t)) + beta[N-1],
01028 CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
01029 else
01030 beta_new[i]=CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
01031 #else //USE_LOGSUMARRAY
01032 float64_t sum=-CMath::INFTY;
01033 for (j=0; j<N; j++)
01034 sum= CMath::logarithmic_sum(sum, get_a(i, j) + get_b(j, p_observations->get_feature(dimension,t)) + beta[j]);
01035
01036 beta_new[i]=sum;
01037 #endif //USE_LOGSUMARRAY
01038 }
01039
01040 if (!BETA_CACHE(dimension).table)
01041 {
01042 dummy=beta;
01043 beta=beta_new;
01044 beta_new=dummy;
01045 }
01046 else
01047 {
01048 beta=beta_new;
01049 beta_new-=N;
01050 }
01051 }
01052
01053 if (time>=0)
01054 {
01055 register int32_t j ;
01056 #ifdef USE_LOGSUMARRAY
01057 for (j=0; j<(N>>1); j++)
01058 ARRAYS(dimension)[j]=CMath::logarithmic_sum(
01059 get_a(state, j<<1) + get_b(j<<1, p_observations->get_feature(dimension,time+1)) + beta[j<<1],
01060 get_a(state, (j<<1)+1) + get_b((j<<1)+1, p_observations->get_feature(dimension,time+1)) + beta[(j<<1)+1]);
01061 if (N%2==1)
01062 return CMath::logarithmic_sum(get_a(state, N-1) + get_b(N-1, p_observations->get_feature(dimension,time+1)) + beta[N-1],
01063 CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
01064 else
01065 return CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
01066 #else //USE_LOGSUMARRAY
01067 float64_t sum=-CMath::INFTY;
01068 for (j=0; j<N; j++)
01069 sum= CMath::logarithmic_sum(sum, get_a(state, j) + get_b(j, p_observations->get_feature(dimension,time+1))+beta[j]);
01070
01071 return sum;
01072 #endif //USE_LOGSUMARRAY
01073 }
01074 else
01075 {
01076 if (BETA_CACHE(dimension).table)
01077 {
01078 #ifdef USE_LOGSUMARRAY//AAA
01079 for (int32_t j=0; j<(N>>1); j++)
01080 ARRAYS(dimension)[j]=CMath::logarithmic_sum(get_p(j<<1) + get_b(j<<1, p_observations->get_feature(dimension,0))+beta[j<<1],
01081 get_p((j<<1)+1) + get_b((j<<1)+1, p_observations->get_feature(dimension,0))+beta[(j<<1)+1]) ;
01082 if (N%2==1)
01083 BETA_CACHE(dimension).sum=CMath::logarithmic_sum(get_p(N-1) + get_b(N-1, p_observations->get_feature(dimension,0))+beta[N-1],
01084 CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
01085 else
01086 BETA_CACHE(dimension).sum=CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
01087 #else //USE_LOGSUMARRAY
01088 float64_t sum=-CMath::INFTY;
01089 for (register int32_t j=0; j<N; j++)
01090 sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
01091 BETA_CACHE(dimension).sum=sum;
01092 #endif //USE_LOGSUMARRAY
01093 BETA_CACHE(dimension).dimension=dimension;
01094 BETA_CACHE(dimension).updated=true;
01095
01096 if (wanted_time<p_observations->get_vector_length(dimension))
01097 return BETA_CACHE(dimension).table[wanted_time*N+state];
01098 else
01099 return BETA_CACHE(dimension).sum;
01100 }
01101 else
01102 {
01103 float64_t sum=-CMath::INFTY;
01104 for (register int32_t j=0; j<N; j++)
01105 sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
01106 return sum;
01107 }
01108 }
01109 }
01110 }
01111
01112
01113
01114 float64_t CHMM::best_path(int32_t dimension)
01115 {
01116 if (!p_observations)
01117 return -1;
01118
01119 if (dimension==-1)
01120 {
01121 if (!all_path_prob_updated)
01122 {
01123 SG_INFO( "computing full viterbi likelihood\n") ;
01124 float64_t sum = 0 ;
01125 for (int32_t i=0; i<p_observations->get_num_vectors(); i++)
01126 sum+=best_path(i) ;
01127 sum /= p_observations->get_num_vectors() ;
01128 all_pat_prob=sum ;
01129 all_path_prob_updated=true ;
01130 return sum ;
01131 } else
01132 return all_pat_prob ;
01133 } ;
01134
01135 if (!STATES_PER_OBSERVATION_PSI(dimension))
01136 return -1 ;
01137
01138 if (dimension >= p_observations->get_num_vectors())
01139 return -1;
01140
01141 if (PATH_PROB_UPDATED(dimension) && dimension==PATH_PROB_DIMENSION(dimension))
01142 return pat_prob;
01143 else
01144 {
01145 register float64_t* delta= ARRAYN2(dimension);
01146 register float64_t* delta_new= ARRAYN1(dimension);
01147
01148 {
01149 for (register int32_t i=0; i<N; i++)
01150 {
01151 delta[i]=get_p(i)+get_b(i, p_observations->get_feature(dimension,0));
01152 set_psi(0, i, 0, dimension);
01153 }
01154 }
01155
01156 #ifdef USE_PATHDEBUG
01157 float64_t worst=-CMath::INFTY/4 ;
01158 #endif
01159
01160 for (register int32_t t=1; t<p_observations->get_vector_length(dimension); t++)
01161 {
01162 register float64_t* dummy;
01163 register int32_t NN=N ;
01164 for (register int32_t j=0; j<NN; j++)
01165 {
01166 register float64_t * matrix_a=&transition_matrix_a[j*N] ;
01167 register float64_t maxj=delta[0] + matrix_a[0];
01168 register int32_t argmax=0;
01169
01170 for (register int32_t i=1; i<NN; i++)
01171 {
01172 register float64_t temp = delta[i] + matrix_a[i];
01173
01174 if (temp>maxj)
01175 {
01176 maxj=temp;
01177 argmax=i;
01178 }
01179 }
01180 #ifdef FIX_POS
01181 if ((!model) || (model->get_fix_pos_state(t,j,NN)!=Model::FIX_DISALLOWED))
01182 #endif
01183 delta_new[j]=maxj + get_b(j,p_observations->get_feature(dimension,t));
01184 #ifdef FIX_POS
01185 else
01186 delta_new[j]=maxj + get_b(j,p_observations->get_feature(dimension,t)) + Model::DISALLOWED_PENALTY;
01187 #endif
01188 set_psi(t, j, argmax, dimension);
01189 }
01190
01191 #ifdef USE_PATHDEBUG
01192 float64_t best=log(0) ;
01193 for (int32_t jj=0; jj<N; jj++)
01194 if (delta_new[jj]>best)
01195 best=delta_new[jj] ;
01196
01197 if (best<-CMath::INFTY/2)
01198 {
01199 SG_DEBUG( "worst case at %i: %e:%e\n", t, best, worst) ;
01200 worst=best ;
01201 } ;
01202 #endif
01203
01204 dummy=delta;
01205 delta=delta_new;
01206 delta_new=dummy;
01207 }
01208
01209
01210 {
01211 register float64_t maxj=delta[0]+get_q(0);
01212 register int32_t argmax=0;
01213
01214 for (register int32_t i=1; i<N; i++)
01215 {
01216 register float64_t temp=delta[i]+get_q(i);
01217
01218 if (temp>maxj)
01219 {
01220 maxj=temp;
01221 argmax=i;
01222 }
01223 }
01224 pat_prob=maxj;
01225 PATH(dimension)[p_observations->get_vector_length(dimension)-1]=argmax;
01226 } ;
01227
01228
01229 {
01230 for (register int32_t t=p_observations->get_vector_length(dimension)-1; t>0; t--)
01231 {
01232 PATH(dimension)[t-1]=get_psi(t, PATH(dimension)[t], dimension);
01233 }
01234 }
01235 PATH_PROB_UPDATED(dimension)=true;
01236 PATH_PROB_DIMENSION(dimension)=dimension;
01237 return pat_prob ;
01238 }
01239 }
01240
01241 #ifndef USE_HMMPARALLEL
01242 float64_t CHMM::model_probability_comp()
01243 {
01244
01245 mod_prob=0 ;
01246 for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
01247 mod_prob+=forward(p_observations->get_vector_length(dim), 0, dim);
01248
01249 mod_prob_updated=true;
01250 return mod_prob;
01251 }
01252
01253 #else
01254
01255 float64_t CHMM::model_probability_comp()
01256 {
01257 pthread_t *threads=SG_MALLOC(pthread_t, parallel->get_num_threads());
01258 S_BW_THREAD_PARAM *params=SG_MALLOC(S_BW_THREAD_PARAM, parallel->get_num_threads());
01259
01260 SG_INFO( "computing full model probablity\n");
01261 mod_prob=0;
01262
01263 for (int32_t cpu=0; cpu<parallel->get_num_threads(); cpu++)
01264 {
01265 params[cpu].hmm=this ;
01266 params[cpu].dim_start= p_observations->get_num_vectors()*cpu/parallel->get_num_threads();
01267 params[cpu].dim_stop= p_observations->get_num_vectors()*(cpu+1)/parallel->get_num_threads();
01268 params[cpu].p_buf=SG_MALLOC(float64_t, N);
01269 params[cpu].q_buf=SG_MALLOC(float64_t, N);
01270 params[cpu].a_buf=SG_MALLOC(float64_t, N*N);
01271 params[cpu].b_buf=SG_MALLOC(float64_t, N*M);
01272 pthread_create(&threads[cpu], NULL, bw_dim_prefetch, (void*)¶ms[cpu]);
01273 }
01274
01275 for (int32_t cpu=0; cpu<parallel->get_num_threads(); cpu++)
01276 {
01277 pthread_join(threads[cpu], NULL);
01278 mod_prob+=params[cpu].ret;
01279 }
01280
01281 for (int32_t i=0; i<parallel->get_num_threads(); i++)
01282 {
01283 SG_FREE(params[i].p_buf);
01284 SG_FREE(params[i].q_buf);
01285 SG_FREE(params[i].a_buf);
01286 SG_FREE(params[i].b_buf);
01287 }
01288
01289 SG_FREE(threads);
01290 SG_FREE(params);
01291
01292 mod_prob_updated=true;
01293 return mod_prob;
01294 }
01295
01296 void* CHMM::bw_dim_prefetch(void* params)
01297 {
01298 CHMM* hmm=((S_BW_THREAD_PARAM*) params)->hmm;
01299 int32_t start=((S_BW_THREAD_PARAM*) params)->dim_start;
01300 int32_t stop=((S_BW_THREAD_PARAM*) params)->dim_stop;
01301 float64_t* p_buf=((S_BW_THREAD_PARAM*) params)->p_buf;
01302 float64_t* q_buf=((S_BW_THREAD_PARAM*) params)->q_buf;
01303 float64_t* a_buf=((S_BW_THREAD_PARAM*) params)->a_buf;
01304 float64_t* b_buf=((S_BW_THREAD_PARAM*) params)->b_buf;
01305 ((S_BW_THREAD_PARAM*)params)->ret=0;
01306
01307 for (int32_t dim=start; dim<stop; dim++)
01308 {
01309 hmm->forward_comp(hmm->p_observations->get_vector_length(dim), hmm->N-1, dim) ;
01310 hmm->backward_comp(hmm->p_observations->get_vector_length(dim), hmm->N-1, dim) ;
01311 float64_t modprob=hmm->model_probability(dim) ;
01312 hmm->ab_buf_comp(p_buf, q_buf, a_buf, b_buf, dim) ;
01313 ((S_BW_THREAD_PARAM*)params)->ret+= modprob;
01314 }
01315 return NULL ;
01316 }
01317
01318 void* CHMM::bw_single_dim_prefetch(void * params)
01319 {
01320 CHMM* hmm=((S_BW_THREAD_PARAM*)params)->hmm ;
01321 int32_t dim=((S_DIM_THREAD_PARAM*)params)->dim ;
01322 ((S_DIM_THREAD_PARAM*)params)->prob_sum = hmm->model_probability(dim);
01323 return NULL ;
01324 }
01325
01326 void* CHMM::vit_dim_prefetch(void * params)
01327 {
01328 CHMM* hmm=((S_DIM_THREAD_PARAM*)params)->hmm ;
01329 int32_t dim=((S_DIM_THREAD_PARAM*)params)->dim ;
01330 ((S_DIM_THREAD_PARAM*)params)->prob_sum = hmm->best_path(dim);
01331 return NULL ;
01332 }
01333
01334 #endif //USE_HMMPARALLEL
01335
01336
01337 #ifdef USE_HMMPARALLEL
01338
01339 void CHMM::ab_buf_comp(
01340 float64_t* p_buf, float64_t* q_buf, float64_t *a_buf, float64_t* b_buf,
01341 int32_t dim)
01342 {
01343 int32_t i,j,t ;
01344 float64_t a_sum;
01345 float64_t b_sum;
01346
01347 float64_t dimmodprob=model_probability(dim);
01348
01349 for (i=0; i<N; i++)
01350 {
01351
01352 p_buf[i]=get_p(i)+get_b(i,p_observations->get_feature(dim,0))+backward(0,i,dim) - dimmodprob;
01353 q_buf[i]=forward(p_observations->get_vector_length(dim)-1, i, dim)+get_q(i) - dimmodprob;
01354
01355
01356 for (j=0; j<N; j++)
01357 {
01358 a_sum=-CMath::INFTY;
01359
01360 for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
01361 {
01362 a_sum= CMath::logarithmic_sum(a_sum, forward(t,i,dim)+
01363 get_a(i,j)+get_b(j,p_observations->get_feature(dim,t+1))+backward(t+1,j,dim));
01364 }
01365 a_buf[N*i+j]=a_sum-dimmodprob ;
01366 }
01367
01368
01369 for (j=0; j<M; j++)
01370 {
01371 b_sum=-CMath::INFTY;
01372
01373 for (t=0; t<p_observations->get_vector_length(dim); t++)
01374 {
01375 if (p_observations->get_feature(dim,t)==j)
01376 b_sum=CMath::logarithmic_sum(b_sum, forward(t,i,dim)+backward(t, i, dim));
01377 }
01378
01379 b_buf[M*i+j]=b_sum-dimmodprob ;
01380 }
01381 }
01382 }
01383
01384
01385 void CHMM::estimate_model_baum_welch(CHMM* hmm)
01386 {
01387 int32_t i,j,cpu;
01388 float64_t fullmodprob=0;
01389
01390
01391 for (i=0; i<N; i++)
01392 {
01393 if (hmm->get_p(i)>CMath::ALMOST_NEG_INFTY)
01394 set_p(i,log(PSEUDO));
01395 else
01396 set_p(i,hmm->get_p(i));
01397 if (hmm->get_q(i)>CMath::ALMOST_NEG_INFTY)
01398 set_q(i,log(PSEUDO));
01399 else
01400 set_q(i,hmm->get_q(i));
01401
01402 for (j=0; j<N; j++)
01403 if (hmm->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
01404 set_a(i,j, log(PSEUDO));
01405 else
01406 set_a(i,j,hmm->get_a(i,j));
01407 for (j=0; j<M; j++)
01408 if (hmm->get_b(i,j)>CMath::ALMOST_NEG_INFTY)
01409 set_b(i,j, log(PSEUDO));
01410 else
01411 set_b(i,j,hmm->get_b(i,j));
01412 }
01413 invalidate_model();
01414
01415 int32_t num_threads = parallel->get_num_threads();
01416
01417 pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
01418 S_BW_THREAD_PARAM *params=SG_MALLOC(S_BW_THREAD_PARAM, num_threads);
01419
01420 if (p_observations->get_num_vectors()<num_threads)
01421 num_threads=p_observations->get_num_vectors();
01422
01423 for (cpu=0; cpu<num_threads; cpu++)
01424 {
01425 params[cpu].p_buf=SG_MALLOC(float64_t, N);
01426 params[cpu].q_buf=SG_MALLOC(float64_t, N);
01427 params[cpu].a_buf=SG_MALLOC(float64_t, N*N);
01428 params[cpu].b_buf=SG_MALLOC(float64_t, N*M);
01429
01430 params[cpu].hmm=hmm;
01431 int32_t start = p_observations->get_num_vectors()*cpu / num_threads;
01432 int32_t stop=p_observations->get_num_vectors()*(cpu+1) / num_threads;
01433
01434 if (cpu == parallel->get_num_threads()-1)
01435 stop=p_observations->get_num_vectors();
01436
01437 ASSERT(start<stop);
01438 params[cpu].dim_start=start;
01439 params[cpu].dim_stop=stop;
01440
01441 pthread_create(&threads[cpu], NULL, bw_dim_prefetch, ¶ms[cpu]);
01442 }
01443
01444 for (cpu=0; cpu<num_threads; cpu++)
01445 {
01446 pthread_join(threads[cpu], NULL);
01447
01448 for (i=0; i<N; i++)
01449 {
01450
01451 set_p(i, CMath::logarithmic_sum(get_p(i), params[cpu].p_buf[i]));
01452 set_q(i, CMath::logarithmic_sum(get_q(i), params[cpu].q_buf[i]));
01453
01454
01455 for (j=0; j<N; j++)
01456 set_a(i,j, CMath::logarithmic_sum(get_a(i,j), params[cpu].a_buf[N*i+j]));
01457
01458
01459 for (j=0; j<M; j++)
01460 set_b(i,j, CMath::logarithmic_sum(get_b(i,j), params[cpu].b_buf[M*i+j]));
01461 }
01462
01463 fullmodprob+=params[cpu].ret;
01464
01465 }
01466
01467 for (cpu=0; cpu<num_threads; cpu++)
01468 {
01469 SG_FREE(params[cpu].p_buf);
01470 SG_FREE(params[cpu].q_buf);
01471 SG_FREE(params[cpu].a_buf);
01472 SG_FREE(params[cpu].b_buf);
01473 }
01474
01475 SG_FREE(threads);
01476 SG_FREE(params);
01477
01478
01479 hmm->mod_prob=fullmodprob;
01480 hmm->mod_prob_updated=true ;
01481
01482
01483 normalize();
01484 invalidate_model();
01485 }
01486
01487 #else // USE_HMMPARALLEL
01488
01489
01490 void CHMM::estimate_model_baum_welch(CHMM* estimate)
01491 {
01492 int32_t i,j,t,dim;
01493 float64_t a_sum, b_sum;
01494 float64_t dimmodprob=0;
01495 float64_t fullmodprob=0;
01496
01497
01498 for (i=0; i<N; i++)
01499 {
01500 if (estimate->get_p(i)>CMath::ALMOST_NEG_INFTY)
01501 set_p(i,log(PSEUDO));
01502 else
01503 set_p(i,estimate->get_p(i));
01504 if (estimate->get_q(i)>CMath::ALMOST_NEG_INFTY)
01505 set_q(i,log(PSEUDO));
01506 else
01507 set_q(i,estimate->get_q(i));
01508
01509 for (j=0; j<N; j++)
01510 if (estimate->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
01511 set_a(i,j, log(PSEUDO));
01512 else
01513 set_a(i,j,estimate->get_a(i,j));
01514 for (j=0; j<M; j++)
01515 if (estimate->get_b(i,j)>CMath::ALMOST_NEG_INFTY)
01516 set_b(i,j, log(PSEUDO));
01517 else
01518 set_b(i,j,estimate->get_b(i,j));
01519 }
01520 invalidate_model();
01521
01522
01523 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
01524 {
01525 dimmodprob=estimate->model_probability(dim);
01526 fullmodprob+=dimmodprob ;
01527
01528 for (i=0; i<N; i++)
01529 {
01530
01531 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));
01532 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 ));
01533
01534 int32_t num = trans_list_backward_cnt[i] ;
01535
01536
01537 for (j=0; j<num; j++)
01538 {
01539 int32_t jj = trans_list_backward[i][j] ;
01540 a_sum=-CMath::INFTY;
01541
01542 for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
01543 {
01544 a_sum= CMath::logarithmic_sum(a_sum, estimate->forward(t,i,dim)+
01545 estimate->get_a(i,jj)+estimate->get_b(jj,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,jj,dim));
01546 }
01547 set_a(i,jj, CMath::logarithmic_sum(get_a(i,jj), a_sum-dimmodprob));
01548 }
01549
01550
01551 for (j=0; j<M; j++)
01552 {
01553 b_sum=-CMath::INFTY;
01554
01555 for (t=0; t<p_observations->get_vector_length(dim); t++)
01556 {
01557 if (p_observations->get_feature(dim,t)==j)
01558 b_sum=CMath::logarithmic_sum(b_sum, estimate->forward(t,i,dim)+estimate->backward(t, i, dim));
01559 }
01560
01561 set_b(i,j, CMath::logarithmic_sum(get_b(i,j), b_sum-dimmodprob));
01562 }
01563 }
01564 }
01565
01566
01567 estimate->mod_prob=fullmodprob;
01568 estimate->mod_prob_updated=true ;
01569
01570
01571 normalize();
01572 invalidate_model();
01573 }
01574
01575
01576 void CHMM::estimate_model_baum_welch_old(CHMM* estimate)
01577 {
01578 int32_t i,j,t,dim;
01579 float64_t a_sum, b_sum;
01580 float64_t dimmodprob=0;
01581 float64_t fullmodprob=0;
01582
01583
01584 for (i=0; i<N; i++)
01585 {
01586 if (estimate->get_p(i)>CMath::ALMOST_NEG_INFTY)
01587 set_p(i,log(PSEUDO));
01588 else
01589 set_p(i,estimate->get_p(i));
01590 if (estimate->get_q(i)>CMath::ALMOST_NEG_INFTY)
01591 set_q(i,log(PSEUDO));
01592 else
01593 set_q(i,estimate->get_q(i));
01594
01595 for (j=0; j<N; j++)
01596 if (estimate->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
01597 set_a(i,j, log(PSEUDO));
01598 else
01599 set_a(i,j,estimate->get_a(i,j));
01600 for (j=0; j<M; j++)
01601 if (estimate->get_b(i,j)>CMath::ALMOST_NEG_INFTY)
01602 set_b(i,j, log(PSEUDO));
01603 else
01604 set_b(i,j,estimate->get_b(i,j));
01605 }
01606 invalidate_model();
01607
01608
01609 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
01610 {
01611 dimmodprob=estimate->model_probability(dim);
01612 fullmodprob+=dimmodprob ;
01613
01614 for (i=0; i<N; i++)
01615 {
01616
01617 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));
01618 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 ));
01619
01620
01621 for (j=0; j<N; j++)
01622 {
01623 a_sum=-CMath::INFTY;
01624
01625 for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
01626 {
01627 a_sum= CMath::logarithmic_sum(a_sum, estimate->forward(t,i,dim)+
01628 estimate->get_a(i,j)+estimate->get_b(j,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,j,dim));
01629 }
01630 set_a(i,j, CMath::logarithmic_sum(get_a(i,j), a_sum-dimmodprob));
01631 }
01632
01633
01634 for (j=0; j<M; j++)
01635 {
01636 b_sum=-CMath::INFTY;
01637
01638 for (t=0; t<p_observations->get_vector_length(dim); t++)
01639 {
01640 if (p_observations->get_feature(dim,t)==j)
01641 b_sum=CMath::logarithmic_sum(b_sum, estimate->forward(t,i,dim)+estimate->backward(t, i, dim));
01642 }
01643
01644 set_b(i,j, CMath::logarithmic_sum(get_b(i,j), b_sum-dimmodprob));
01645 }
01646 }
01647 }
01648
01649
01650 estimate->mod_prob=fullmodprob;
01651 estimate->mod_prob_updated=true ;
01652
01653
01654 normalize();
01655 invalidate_model();
01656 }
01657 #endif // USE_HMMPARALLEL
01658
01659
01660
01661 void CHMM::estimate_model_baum_welch_trans(CHMM* estimate)
01662 {
01663 int32_t i,j,t,dim;
01664 float64_t a_sum;
01665 float64_t dimmodprob=0;
01666 float64_t fullmodprob=0;
01667
01668
01669 for (i=0; i<N; i++)
01670 {
01671 if (estimate->get_p(i)>CMath::ALMOST_NEG_INFTY)
01672 set_p(i,log(PSEUDO));
01673 else
01674 set_p(i,estimate->get_p(i));
01675 if (estimate->get_q(i)>CMath::ALMOST_NEG_INFTY)
01676 set_q(i,log(PSEUDO));
01677 else
01678 set_q(i,estimate->get_q(i));
01679
01680 for (j=0; j<N; j++)
01681 if (estimate->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
01682 set_a(i,j, log(PSEUDO));
01683 else
01684 set_a(i,j,estimate->get_a(i,j));
01685 for (j=0; j<M; j++)
01686 set_b(i,j,estimate->get_b(i,j));
01687 }
01688 invalidate_model();
01689
01690
01691 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
01692 {
01693 dimmodprob=estimate->model_probability(dim);
01694 fullmodprob+=dimmodprob ;
01695
01696 for (i=0; i<N; i++)
01697 {
01698
01699 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));
01700 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 ));
01701
01702 int32_t num = trans_list_backward_cnt[i] ;
01703
01704 for (j=0; j<num; j++)
01705 {
01706 int32_t jj = trans_list_backward[i][j] ;
01707 a_sum=-CMath::INFTY;
01708
01709 for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
01710 {
01711 a_sum= CMath::logarithmic_sum(a_sum, estimate->forward(t,i,dim)+
01712 estimate->get_a(i,jj)+estimate->get_b(jj,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,jj,dim));
01713 }
01714 set_a(i,jj, CMath::logarithmic_sum(get_a(i,jj), a_sum-dimmodprob));
01715 }
01716 }
01717 }
01718
01719
01720 estimate->mod_prob=fullmodprob;
01721 estimate->mod_prob_updated=true ;
01722
01723
01724 normalize();
01725 invalidate_model();
01726 }
01727
01728
01729
01730
01731 void CHMM::estimate_model_baum_welch_defined(CHMM* estimate)
01732 {
01733 int32_t i,j,old_i,k,t,dim;
01734 float64_t a_sum_num, b_sum_num;
01735 float64_t a_sum_denom, b_sum_denom;
01736 float64_t dimmodprob=-CMath::INFTY;
01737 float64_t fullmodprob=0;
01738 float64_t* A=ARRAYN1(0);
01739 float64_t* B=ARRAYN2(0);
01740
01741
01742
01743 for (k=0; (i=model->get_learn_p(k))!=-1; k++)
01744 set_p(i,log(PSEUDO));
01745
01746 for (k=0; (i=model->get_learn_q(k))!=-1; k++)
01747 set_q(i,log(PSEUDO));
01748
01749 for (k=0; (i=model->get_learn_a(k,0))!=-1; k++)
01750 {
01751 j=model->get_learn_a(k,1);
01752 set_a(i,j, log(PSEUDO));
01753 }
01754
01755 for (k=0; (i=model->get_learn_b(k,0))!=-1; k++)
01756 {
01757 j=model->get_learn_b(k,1);
01758 set_b(i,j, log(PSEUDO));
01759 }
01760
01761 for (i=0; i<N; i++)
01762 {
01763 A[i]=log(PSEUDO);
01764 B[i]=log(PSEUDO);
01765 }
01766
01767 #ifdef USE_HMMPARALLEL
01768 int32_t num_threads = parallel->get_num_threads();
01769 pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
01770 S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);
01771
01772 if (p_observations->get_num_vectors()<num_threads)
01773 num_threads=p_observations->get_num_vectors();
01774 #endif
01775
01776
01777 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
01778 {
01779 #ifdef USE_HMMPARALLEL
01780 if (dim%num_threads==0)
01781 {
01782 for (i=0; i<num_threads; i++)
01783 {
01784 if (dim+i<p_observations->get_num_vectors())
01785 {
01786 params[i].hmm=estimate ;
01787 params[i].dim=dim+i ;
01788 pthread_create(&threads[i], NULL, bw_single_dim_prefetch, (void*)¶ms[i]) ;
01789 }
01790 }
01791 for (i=0; i<num_threads; i++)
01792 {
01793 if (dim+i<p_observations->get_num_vectors())
01794 {
01795 pthread_join(threads[i], NULL);
01796 dimmodprob = params[i].prob_sum;
01797 }
01798 }
01799 }
01800 #else
01801 dimmodprob=estimate->model_probability(dim);
01802 #endif // USE_HMMPARALLEL
01803
01804
01805 fullmodprob+= dimmodprob;
01806
01807
01808 for (k=0; (i=model->get_learn_p(k))!=-1; k++)
01809 set_p(i, CMath::logarithmic_sum(get_p(i), estimate->forward(0,i,dim)+estimate->backward(0,i,dim) - dimmodprob ) );
01810
01811 for (k=0; (i=model->get_learn_q(k))!=-1; k++)
01812 set_q(i, CMath::logarithmic_sum(get_q(i), estimate->forward(p_observations->get_vector_length(dim)-1, i, dim)+
01813 estimate->backward(p_observations->get_vector_length(dim)-1, i, dim) - dimmodprob ) );
01814
01815
01816 old_i=-1;
01817 for (k=0; (i=model->get_learn_a(k,0))!=-1; k++)
01818 {
01819
01820
01821 if (old_i!=i)
01822 {
01823 old_i=i;
01824 a_sum_denom=-CMath::INFTY;
01825
01826 for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
01827 a_sum_denom= CMath::logarithmic_sum(a_sum_denom, estimate->forward(t,i,dim)+estimate->backward(t,i,dim));
01828
01829 A[i]= CMath::logarithmic_sum(A[i], a_sum_denom-dimmodprob);
01830 }
01831
01832 j=model->get_learn_a(k,1);
01833 a_sum_num=-CMath::INFTY;
01834 for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
01835 {
01836 a_sum_num= CMath::logarithmic_sum(a_sum_num, estimate->forward(t,i,dim)+
01837 estimate->get_a(i,j)+estimate->get_b(j,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,j,dim));
01838 }
01839
01840 set_a(i,j, CMath::logarithmic_sum(get_a(i,j), a_sum_num-dimmodprob));
01841 }
01842
01843
01844 old_i=-1;
01845 for (k=0; (i=model->get_learn_b(k,0))!=-1; k++)
01846 {
01847
01848
01849
01850 if (old_i!=i)
01851 {
01852 old_i=i;
01853 b_sum_denom=-CMath::INFTY;
01854
01855 for (t=0; t<p_observations->get_vector_length(dim); t++)
01856 b_sum_denom= CMath::logarithmic_sum(b_sum_denom, estimate->forward(t,i,dim)+estimate->backward(t,i,dim));
01857
01858 B[i]= CMath::logarithmic_sum(B[i], b_sum_denom-dimmodprob);
01859 }
01860
01861 j=model->get_learn_b(k,1);
01862 b_sum_num=-CMath::INFTY;
01863 for (t=0; t<p_observations->get_vector_length(dim); t++)
01864 {
01865 if (p_observations->get_feature(dim,t)==j)
01866 b_sum_num=CMath::logarithmic_sum(b_sum_num, estimate->forward(t,i,dim)+estimate->backward(t, i, dim));
01867 }
01868
01869 set_b(i,j, CMath::logarithmic_sum(get_b(i,j), b_sum_num-dimmodprob));
01870 }
01871 }
01872 #ifdef USE_HMMPARALLEL
01873 SG_FREE(threads);
01874 SG_FREE(params);
01875 #endif
01876
01877
01878
01879 for (k=0; (i=model->get_learn_p(k))!=-1; k++)
01880 set_p(i, get_p(i)-log(p_observations->get_num_vectors()+N*PSEUDO) );
01881
01882 for (k=0; (i=model->get_learn_q(k))!=-1; k++)
01883 set_q(i, get_q(i)-log(p_observations->get_num_vectors()+N*PSEUDO) );
01884
01885 for (k=0; (i=model->get_learn_a(k,0))!=-1; k++)
01886 {
01887 j=model->get_learn_a(k,1);
01888 set_a(i,j, get_a(i,j) - A[i]);
01889 }
01890
01891 for (k=0; (i=model->get_learn_b(k,0))!=-1; k++)
01892 {
01893 j=model->get_learn_b(k,1);
01894 set_b(i,j, get_b(i,j) - B[i]);
01895 }
01896
01897
01898 estimate->mod_prob=fullmodprob;
01899 estimate->mod_prob_updated=true ;
01900
01901
01902 normalize();
01903 invalidate_model();
01904 }
01905
01906
01907 void CHMM::estimate_model_viterbi(CHMM* estimate)
01908 {
01909 int32_t i,j,t;
01910 float64_t sum;
01911 float64_t* P=ARRAYN1(0);
01912 float64_t* Q=ARRAYN2(0);
01913
01914 path_deriv_updated=false ;
01915
01916
01917 for (i=0; i<N; i++)
01918 {
01919 for (j=0; j<N; j++)
01920 set_A(i,j, PSEUDO);
01921
01922 for (j=0; j<M; j++)
01923 set_B(i,j, PSEUDO);
01924
01925 P[i]=PSEUDO;
01926 Q[i]=PSEUDO;
01927 }
01928
01929 float64_t allpatprob=0 ;
01930
01931 #ifdef USE_HMMPARALLEL
01932 int32_t num_threads = parallel->get_num_threads();
01933 pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
01934 S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);
01935
01936 if (p_observations->get_num_vectors()<num_threads)
01937 num_threads=p_observations->get_num_vectors();
01938 #endif
01939
01940 for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
01941 {
01942
01943 #ifdef USE_HMMPARALLEL
01944 if (dim%num_threads==0)
01945 {
01946 for (i=0; i<num_threads; i++)
01947 {
01948 if (dim+i<p_observations->get_num_vectors())
01949 {
01950 params[i].hmm=estimate ;
01951 params[i].dim=dim+i ;
01952 pthread_create(&threads[i], NULL, vit_dim_prefetch, (void*)¶ms[i]) ;
01953 }
01954 }
01955 for (i=0; i<num_threads; i++)
01956 {
01957 if (dim+i<p_observations->get_num_vectors())
01958 {
01959 pthread_join(threads[i], NULL);
01960 allpatprob += params[i].prob_sum;
01961 }
01962 }
01963 }
01964 #else
01965
01966 allpatprob += estimate->best_path(dim);
01967 #endif // USE_HMMPARALLEL
01968
01969
01970 for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
01971 {
01972 set_A(estimate->PATH(dim)[t], estimate->PATH(dim)[t+1], get_A(estimate->PATH(dim)[t], estimate->PATH(dim)[t+1])+1);
01973 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);
01974 }
01975
01976 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 );
01977
01978 P[estimate->PATH(dim)[0]]++;
01979 Q[estimate->PATH(dim)[p_observations->get_vector_length(dim)-1]]++;
01980 }
01981
01982 #ifdef USE_HMMPARALLEL
01983 SG_FREE(threads);
01984 SG_FREE(params);
01985 #endif
01986
01987 allpatprob/=p_observations->get_num_vectors() ;
01988 estimate->all_pat_prob=allpatprob ;
01989 estimate->all_path_prob_updated=true ;
01990
01991
01992 for (i=0; i<N; i++)
01993 {
01994 sum=0;
01995 for (j=0; j<N; j++)
01996 sum+=get_A(i,j);
01997
01998 for (j=0; j<N; j++)
01999 set_a(i,j, log(get_A(i,j)/sum));
02000 }
02001
02002
02003 for (i=0; i<N; i++)
02004 {
02005 sum=0;
02006 for (j=0; j<M; j++)
02007 sum+=get_B(i,j);
02008
02009 for (j=0; j<M; j++)
02010 set_b(i,j, log(get_B(i, j)/sum));
02011 }
02012
02013
02014 sum=0;
02015 for (i=0; i<N; i++)
02016 sum+=P[i];
02017
02018 for (i=0; i<N; i++)
02019 set_p(i, log(P[i]/sum));
02020
02021
02022 sum=0;
02023 for (i=0; i<N; i++)
02024 sum+=Q[i];
02025
02026 for (i=0; i<N; i++)
02027 set_q(i, log(Q[i]/sum));
02028
02029
02030 invalidate_model();
02031 }
02032
02033
02034 void CHMM::estimate_model_viterbi_defined(CHMM* estimate)
02035 {
02036 int32_t i,j,k,t;
02037 float64_t sum;
02038 float64_t* P=ARRAYN1(0);
02039 float64_t* Q=ARRAYN2(0);
02040
02041 path_deriv_updated=false ;
02042
02043
02044 for (i=0; i<N; i++)
02045 {
02046 for (j=0; j<N; j++)
02047 set_A(i,j, PSEUDO);
02048
02049 for (j=0; j<M; j++)
02050 set_B(i,j, PSEUDO);
02051
02052 P[i]=PSEUDO;
02053 Q[i]=PSEUDO;
02054 }
02055
02056 #ifdef USE_HMMPARALLEL
02057 int32_t num_threads = parallel->get_num_threads();
02058 pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
02059 S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);
02060 #endif
02061
02062 float64_t allpatprob=0.0 ;
02063 for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
02064 {
02065
02066 #ifdef USE_HMMPARALLEL
02067 if (dim%num_threads==0)
02068 {
02069 for (i=0; i<num_threads; i++)
02070 {
02071 if (dim+i<p_observations->get_num_vectors())
02072 {
02073 params[i].hmm=estimate ;
02074 params[i].dim=dim+i ;
02075 pthread_create(&threads[i], NULL, vit_dim_prefetch, (void*)¶ms[i]) ;
02076 }
02077 }
02078 for (i=0; i<num_threads; i++)
02079 {
02080 if (dim+i<p_observations->get_num_vectors())
02081 {
02082 pthread_join(threads[i], NULL);
02083 allpatprob += params[i].prob_sum;
02084 }
02085 }
02086 }
02087 #else // USE_HMMPARALLEL
02088
02089 allpatprob += estimate->best_path(dim);
02090 #endif // USE_HMMPARALLEL
02091
02092
02093
02094 for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
02095 {
02096 set_A(estimate->PATH(dim)[t], estimate->PATH(dim)[t+1], get_A(estimate->PATH(dim)[t], estimate->PATH(dim)[t+1])+1);
02097 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);
02098 }
02099
02100 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 );
02101
02102 P[estimate->PATH(dim)[0]]++;
02103 Q[estimate->PATH(dim)[p_observations->get_vector_length(dim)-1]]++;
02104 }
02105
02106 #ifdef USE_HMMPARALLEL
02107 SG_FREE(threads);
02108 SG_FREE(params);
02109 #endif
02110
02111
02112
02113
02114 allpatprob/=p_observations->get_num_vectors() ;
02115 estimate->all_pat_prob=allpatprob ;
02116 estimate->all_path_prob_updated=true ;
02117
02118
02119
02120 for (i=0; i<N; i++)
02121 {
02122 for (j=0; j<N; j++)
02123 set_a(i,j, estimate->get_a(i,j));
02124
02125 for (j=0; j<M; j++)
02126 set_b(i,j, estimate->get_b(i,j));
02127 }
02128
02129
02130 i=0;
02131 sum=0;
02132 j=model->get_learn_a(i,0);
02133 k=i;
02134 while (model->get_learn_a(i,0)!=-1 || k<i)
02135 {
02136 if (j==model->get_learn_a(i,0))
02137 {
02138 sum+=get_A(model->get_learn_a(i,0), model->get_learn_a(i,1));
02139 i++;
02140 }
02141 else
02142 {
02143 while (k<i)
02144 {
02145 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));
02146 k++;
02147 }
02148
02149 sum=0;
02150 j=model->get_learn_a(i,0);
02151 k=i;
02152 }
02153 }
02154
02155
02156 i=0;
02157 sum=0;
02158 j=model->get_learn_b(i,0);
02159 k=i;
02160 while (model->get_learn_b(i,0)!=-1 || k<i)
02161 {
02162 if (j==model->get_learn_b(i,0))
02163 {
02164 sum+=get_B(model->get_learn_b(i,0),model->get_learn_b(i,1));
02165 i++;
02166 }
02167 else
02168 {
02169 while (k<i)
02170 {
02171 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));
02172 k++;
02173 }
02174
02175 sum=0;
02176 j=model->get_learn_b(i,0);
02177 k=i;
02178 }
02179 }
02180
02181 i=0;
02182 sum=0;
02183 while (model->get_learn_p(i)!=-1)
02184 {
02185 sum+=P[model->get_learn_p(i)] ;
02186 i++ ;
02187 } ;
02188 i=0 ;
02189 while (model->get_learn_p(i)!=-1)
02190 {
02191 set_p(model->get_learn_p(i), log(P[model->get_learn_p(i)]/sum));
02192 i++ ;
02193 } ;
02194
02195 i=0;
02196 sum=0;
02197 while (model->get_learn_q(i)!=-1)
02198 {
02199 sum+=Q[model->get_learn_q(i)] ;
02200 i++ ;
02201 } ;
02202 i=0 ;
02203 while (model->get_learn_q(i)!=-1)
02204 {
02205 set_q(model->get_learn_q(i), log(Q[model->get_learn_q(i)]/sum));
02206 i++ ;
02207 } ;
02208
02209
02210
02211 invalidate_model();
02212 }
02213
02214
02215
02216 void CHMM::output_model(bool verbose)
02217 {
02218 int32_t i,j;
02219 float64_t checksum;
02220
02221
02222 SG_INFO( "log(Pr[O|model])=%e, #states: %i, #observationssymbols: %i, #observations: %ix%i\n",
02223 (float64_t)((p_observations) ? model_probability() : -CMath::INFTY),
02224 N, M, ((p_observations) ? p_observations->get_max_vector_length() : 0), ((p_observations) ? p_observations->get_num_vectors() : 0));
02225
02226 if (verbose)
02227 {
02228
02229 SG_INFO( "\ntransition matrix\n");
02230 for (i=0; i<N; i++)
02231 {
02232 checksum= get_q(i);
02233 for (j=0; j<N; j++)
02234 {
02235 checksum= CMath::logarithmic_sum(checksum, get_a(i,j));
02236
02237 SG_INFO( "a(%02i,%02i)=%1.4f ",i,j, (float32_t) exp(get_a(i,j)));
02238
02239 if (j % 4 == 3)
02240 SG_PRINT( "\n");
02241 }
02242 if (fabs(checksum)>1e-5)
02243 SG_DEBUG( " checksum % E ******* \n",checksum);
02244 else
02245 SG_DEBUG( " checksum % E\n",checksum);
02246 }
02247
02248
02249 SG_INFO( "\ndistribution of start states\n");
02250 checksum=-CMath::INFTY;
02251 for (i=0; i<N; i++)
02252 {
02253 checksum= CMath::logarithmic_sum(checksum, get_p(i));
02254 SG_INFO( "p(%02i)=%1.4f ",i, (float32_t) exp(get_p(i)));
02255 if (i % 4 == 3)
02256 SG_PRINT( "\n");
02257 }
02258 if (fabs(checksum)>1e-5)
02259 SG_DEBUG( " checksum % E ******* \n",checksum);
02260 else
02261 SG_DEBUG( " checksum=% E\n", checksum);
02262
02263
02264 SG_INFO( "\ndistribution of terminal states\n");
02265 checksum=-CMath::INFTY;
02266 for (i=0; i<N; i++)
02267 {
02268 checksum= CMath::logarithmic_sum(checksum, get_q(i));
02269 SG_INFO("q(%02i)=%1.4f ",i, (float32_t) exp(get_q(i)));
02270 if (i % 4 == 3)
02271 SG_INFO("\n");
02272 }
02273 if (fabs(checksum)>1e-5)
02274 SG_DEBUG( " checksum % E ******* \n",checksum);
02275 else
02276 SG_DEBUG( " checksum=% E\n", checksum);
02277
02278
02279 SG_INFO("\ndistribution of observations given the state\n");
02280 for (i=0; i<N; i++)
02281 {
02282 checksum=-CMath::INFTY;
02283 for (j=0; j<M; j++)
02284 {
02285 checksum=CMath::logarithmic_sum(checksum, get_b(i,j));
02286 SG_INFO("b(%02i,%02i)=%1.4f ",i,j, (float32_t) exp(get_b(i,j)));
02287 if (j % 4 == 3)
02288 SG_PRINT("\n");
02289 }
02290 if (fabs(checksum)>1e-5)
02291 SG_DEBUG( " checksum % E ******* \n",checksum);
02292 else
02293 SG_DEBUG( " checksum % E\n",checksum);
02294 }
02295 }
02296 SG_PRINT("\n");
02297 }
02298
02299
02300 void CHMM::output_model_defined(bool verbose)
02301 {
02302 int32_t i,j;
02303 if (!model)
02304 return ;
02305
02306
02307 SG_INFO("log(Pr[O|model])=%e, #states: %i, #observationssymbols: %i, #observations: %ix%i\n",
02308 (float64_t)((p_observations) ? model_probability() : -CMath::INFTY),
02309 N, M, ((p_observations) ? p_observations->get_max_vector_length() : 0), ((p_observations) ? p_observations->get_num_vectors() : 0));
02310
02311 if (verbose)
02312 {
02313
02314 SG_INFO("\ntransition matrix\n");
02315
02316
02317 i=0;
02318 j=model->get_learn_a(i,0);
02319 while (model->get_learn_a(i,0)!=-1)
02320 {
02321 if (j!=model->get_learn_a(i,0))
02322 {
02323 j=model->get_learn_a(i,0);
02324 SG_PRINT("\n");
02325 }
02326
02327 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))));
02328 i++;
02329 }
02330
02331
02332 SG_INFO("\n\ndistribution of observations given the state\n");
02333 i=0;
02334 j=model->get_learn_b(i,0);
02335 while (model->get_learn_b(i,0)!=-1)
02336 {
02337 if (j!=model->get_learn_b(i,0))
02338 {
02339 j=model->get_learn_b(i,0);
02340 SG_PRINT("\n");
02341 }
02342
02343 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))));
02344 i++;
02345 }
02346
02347 SG_PRINT("\n");
02348 }
02349 SG_PRINT("\n");
02350 }
02351
02352
02353
02354
02355 void CHMM::convert_to_log()
02356 {
02357 int32_t i,j;
02358
02359 for (i=0; i<N; i++)
02360 {
02361 if (get_p(i)!=0)
02362 set_p(i, log(get_p(i)));
02363 else
02364 set_p(i, -CMath::INFTY);;
02365 }
02366
02367 for (i=0; i<N; i++)
02368 {
02369 if (get_q(i)!=0)
02370 set_q(i, log(get_q(i)));
02371 else
02372 set_q(i, -CMath::INFTY);;
02373 }
02374
02375 for (i=0; i<N; i++)
02376 {
02377 for (j=0; j<N; j++)
02378 {
02379 if (get_a(i,j)!=0)
02380 set_a(i,j, log(get_a(i,j)));
02381 else
02382 set_a(i,j, -CMath::INFTY);
02383 }
02384 }
02385
02386 for (i=0; i<N; i++)
02387 {
02388 for (j=0; j<M; j++)
02389 {
02390 if (get_b(i,j)!=0)
02391 set_b(i,j, log(get_b(i,j)));
02392 else
02393 set_b(i,j, -CMath::INFTY);
02394 }
02395 }
02396 loglikelihood=true;
02397
02398 invalidate_model();
02399 }
02400
02401
02402 void CHMM::init_model_random()
02403 {
02404 const float64_t MIN_RAND=23e-3;
02405
02406 float64_t sum;
02407 int32_t i,j;
02408
02409
02410 for (i=0; i<N; i++)
02411 {
02412 sum=0;
02413 for (j=0; j<N; j++)
02414 {
02415 set_a(i,j, CMath::random(MIN_RAND, 1.0));
02416
02417 sum+=get_a(i,j);
02418 }
02419
02420 for (j=0; j<N; j++)
02421 set_a(i,j, get_a(i,j)/sum);
02422 }
02423
02424
02425 sum=0;
02426 for (i=0; i<N; i++)
02427 {
02428 set_p(i, CMath::random(MIN_RAND, 1.0));
02429
02430 sum+=get_p(i);
02431 }
02432
02433 for (i=0; i<N; i++)
02434 set_p(i, get_p(i)/sum);
02435
02436
02437 sum=0;
02438 for (i=0; i<N; i++)
02439 {
02440 set_q(i, CMath::random(MIN_RAND, 1.0));
02441
02442 sum+=get_q(i);
02443 }
02444
02445 for (i=0; i<N; i++)
02446 set_q(i, get_q(i)/sum);
02447
02448
02449 for (i=0; i<N; i++)
02450 {
02451 sum=0;
02452 for (j=0; j<M; j++)
02453 {
02454 set_b(i,j, CMath::random(MIN_RAND, 1.0));
02455
02456 sum+=get_b(i,j);
02457 }
02458
02459 for (j=0; j<M; j++)
02460 set_b(i,j, get_b(i,j)/sum);
02461 }
02462
02463
02464 invalidate_model();
02465 }
02466
02467
02468 void CHMM::init_model_defined()
02469 {
02470 int32_t i,j,k,r;
02471 float64_t sum;
02472 const float64_t MIN_RAND=23e-3;
02473
02474
02475 for (i=0; i<N; i++)
02476 for (j=0; j<N; j++)
02477 set_a(i,j, 0);
02478
02479
02480 for (i=0; i<N; i++)
02481 set_p(i, 0);
02482
02483
02484 for (i=0; i<N; i++)
02485 set_q(i, 0);
02486
02487
02488 for (i=0; i<N; i++)
02489 for (j=0; j<M; j++)
02490 set_b(i,j, 0);
02491
02492
02493
02494 float64_t *R=SG_MALLOC(float64_t, N);
02495 for (r=0; r<N; r++) R[r]=CMath::random(MIN_RAND,1.0);
02496 i=0; sum=0; k=i;
02497 j=model->get_learn_a(i,0);
02498 while (model->get_learn_a(i,0)!=-1 || k<i)
02499 {
02500 if (j==model->get_learn_a(i,0))
02501 {
02502 sum+=R[model->get_learn_a(i,1)] ;
02503 i++;
02504 }
02505 else
02506 {
02507 while (k<i)
02508 {
02509 set_a(model->get_learn_a(k,0), model->get_learn_a(k,1),
02510 R[model->get_learn_a(k,1)]/sum);
02511 k++;
02512 }
02513 j=model->get_learn_a(i,0);
02514 k=i;
02515 sum=0;
02516 for (r=0; r<N; r++) R[r]=CMath::random(MIN_RAND,1.0);
02517 }
02518 }
02519 SG_FREE(R); R=NULL ;
02520
02521
02522 R=SG_MALLOC(float64_t, M);
02523 for (r=0; r<M; r++) R[r]=CMath::random(MIN_RAND,1.0);
02524 i=0; sum=0; k=0 ;
02525 j=model->get_learn_b(i,0);
02526 while (model->get_learn_b(i,0)!=-1 || k<i)
02527 {
02528 if (j==model->get_learn_b(i,0))
02529 {
02530 sum+=R[model->get_learn_b(i,1)] ;
02531 i++;
02532 }
02533 else
02534 {
02535 while (k<i)
02536 {
02537 set_b(model->get_learn_b(k,0),model->get_learn_b(k,1),
02538 R[model->get_learn_b(k,1)]/sum);
02539 k++;
02540 }
02541
02542 j=model->get_learn_b(i,0);
02543 k=i;
02544 sum=0;
02545 for (r=0; r<M; r++) R[r]=CMath::random(MIN_RAND,1.0);
02546 }
02547 }
02548 SG_FREE(R); R=NULL ;
02549
02550
02551 i=0;
02552 while (model->get_const_a(i,0) != -1)
02553 {
02554 set_a(model->get_const_a(i,0), model->get_const_a(i,1), model->get_const_a_val(i));
02555 i++;
02556 }
02557
02558
02559 i=0;
02560 while (model->get_const_b(i,0) != -1)
02561 {
02562 set_b(model->get_const_b(i,0), model->get_const_b(i,1), model->get_const_b_val(i));
02563 i++;
02564 }
02565
02566
02567 i=0;
02568 while (model->get_const_p(i) != -1)
02569 {
02570 set_p(model->get_const_p(i), model->get_const_p_val(i));
02571 i++;
02572 }
02573
02574
02575 for (i=0; i<N; i++)
02576 set_q(i, 0.0);
02577
02578
02579 i=0;
02580 while (model->get_const_q(i) != -1)
02581 {
02582 set_q(model->get_const_q(i), model->get_const_q_val(i));
02583 i++;
02584 }
02585
02586
02587 i=0;
02588 sum=0;
02589 while (model->get_learn_p(i)!=-1)
02590 {
02591 set_p(model->get_learn_p(i),CMath::random(MIN_RAND,1.0)) ;
02592 sum+=get_p(model->get_learn_p(i)) ;
02593 i++ ;
02594 } ;
02595 i=0 ;
02596 while (model->get_learn_p(i)!=-1)
02597 {
02598 set_p(model->get_learn_p(i), get_p(model->get_learn_p(i))/sum);
02599 i++ ;
02600 } ;
02601
02602
02603 i=0;
02604 sum=0;
02605 while (model->get_learn_q(i)!=-1)
02606 {
02607 set_q(model->get_learn_q(i),CMath::random(MIN_RAND,1.0)) ;
02608 sum+=get_q(model->get_learn_q(i)) ;
02609 i++ ;
02610 } ;
02611 i=0 ;
02612 while (model->get_learn_q(i)!=-1)
02613 {
02614 set_q(model->get_learn_q(i), get_q(model->get_learn_q(i))/sum);
02615 i++ ;
02616 } ;
02617
02618
02619 invalidate_model();
02620 }
02621
02622 void CHMM::clear_model()
02623 {
02624 int32_t i,j;
02625 for (i=0; i<N; i++)
02626 {
02627 set_p(i, log(PSEUDO));
02628 set_q(i, log(PSEUDO));
02629
02630 for (j=0; j<N; j++)
02631 set_a(i,j, log(PSEUDO));
02632
02633 for (j=0; j<M; j++)
02634 set_b(i,j, log(PSEUDO));
02635 }
02636 }
02637
02638 void CHMM::clear_model_defined()
02639 {
02640 int32_t i,j,k;
02641
02642 for (i=0; (j=model->get_learn_p(i))!=-1; i++)
02643 set_p(j, log(PSEUDO));
02644
02645 for (i=0; (j=model->get_learn_q(i))!=-1; i++)
02646 set_q(j, log(PSEUDO));
02647
02648 for (i=0; (j=model->get_learn_a(i,0))!=-1; i++)
02649 {
02650 k=model->get_learn_a(i,1);
02651 set_a(j,k, log(PSEUDO));
02652 }
02653
02654 for (i=0; (j=model->get_learn_b(i,0))!=-1; i++)
02655 {
02656 k=model->get_learn_b(i,1);
02657 set_b(j,k, log(PSEUDO));
02658 }
02659 }
02660
02661 void CHMM::copy_model(CHMM* l)
02662 {
02663 int32_t i,j;
02664 for (i=0; i<N; i++)
02665 {
02666 set_p(i, l->get_p(i));
02667 set_q(i, l->get_q(i));
02668
02669 for (j=0; j<N; j++)
02670 set_a(i,j, l->get_a(i,j));
02671
02672 for (j=0; j<M; j++)
02673 set_b(i,j, l->get_b(i,j));
02674 }
02675 }
02676
02677 void CHMM::invalidate_model()
02678 {
02679
02680 this->mod_prob=0.0;
02681 this->mod_prob_updated=false;
02682
02683 if (mem_initialized)
02684 {
02685 if (trans_list_forward_cnt)
02686 SG_FREE(trans_list_forward_cnt);
02687 trans_list_forward_cnt=NULL ;
02688 if (trans_list_backward_cnt)
02689 SG_FREE(trans_list_backward_cnt);
02690 trans_list_backward_cnt=NULL ;
02691 if (trans_list_forward)
02692 {
02693 for (int32_t i=0; i<trans_list_len; i++)
02694 if (trans_list_forward[i])
02695 SG_FREE(trans_list_forward[i]);
02696 SG_FREE(trans_list_forward);
02697 trans_list_forward=NULL ;
02698 }
02699 if (trans_list_backward)
02700 {
02701 for (int32_t i=0; i<trans_list_len; i++)
02702 if (trans_list_backward[i])
02703 SG_FREE(trans_list_backward[i]);
02704 SG_FREE(trans_list_backward);
02705 trans_list_backward = NULL ;
02706 } ;
02707
02708 trans_list_len = N ;
02709 trans_list_forward = SG_MALLOC(T_STATES*, N);
02710 trans_list_forward_cnt = SG_MALLOC(T_STATES, N);
02711
02712 for (int32_t j=0; j<N; j++)
02713 {
02714 trans_list_forward_cnt[j]= 0 ;
02715 trans_list_forward[j]= SG_MALLOC(T_STATES, N);
02716 for (int32_t i=0; i<N; i++)
02717 if (get_a(i,j)>CMath::ALMOST_NEG_INFTY)
02718 {
02719 trans_list_forward[j][trans_list_forward_cnt[j]]=i ;
02720 trans_list_forward_cnt[j]++ ;
02721 }
02722 } ;
02723
02724 trans_list_backward = SG_MALLOC(T_STATES*, N);
02725 trans_list_backward_cnt = SG_MALLOC(T_STATES, N);
02726
02727 for (int32_t i=0; i<N; i++)
02728 {
02729 trans_list_backward_cnt[i]= 0 ;
02730 trans_list_backward[i]= SG_MALLOC(T_STATES, N);
02731 for (int32_t j=0; j<N; j++)
02732 if (get_a(i,j)>CMath::ALMOST_NEG_INFTY)
02733 {
02734 trans_list_backward[i][trans_list_backward_cnt[i]]=j ;
02735 trans_list_backward_cnt[i]++ ;
02736 }
02737 } ;
02738 } ;
02739 this->all_pat_prob=0.0;
02740 this->pat_prob=0.0;
02741 this->path_deriv_updated=false ;
02742 this->path_deriv_dimension=-1 ;
02743 this->all_path_prob_updated=false;
02744
02745 #ifdef USE_HMMPARALLEL_STRUCTURES
02746 {
02747 for (int32_t i=0; i<parallel->get_num_threads(); i++)
02748 {
02749 this->alpha_cache[i].updated=false;
02750 this->beta_cache[i].updated=false;
02751 path_prob_updated[i]=false ;
02752 path_prob_dimension[i]=-1 ;
02753 } ;
02754 }
02755 #else // USE_HMMPARALLEL_STRUCTURES
02756 this->alpha_cache.updated=false;
02757 this->beta_cache.updated=false;
02758 this->path_prob_dimension=-1;
02759 this->path_prob_updated=false;
02760
02761 #endif // USE_HMMPARALLEL_STRUCTURES
02762 }
02763
02764 void CHMM::open_bracket(FILE* file)
02765 {
02766 int32_t value;
02767 while (((value=fgetc(file)) != EOF) && (value!='['))
02768 {
02769 if (value=='\n')
02770 line++;
02771 }
02772
02773 if (value==EOF)
02774 error(line, "expected \"[\" in input file");
02775
02776 while (((value=fgetc(file)) != EOF) && (isspace(value)))
02777 {
02778 if (value=='\n')
02779 line++;
02780 }
02781
02782 ungetc(value, file);
02783 }
02784
02785 void CHMM::close_bracket(FILE* file)
02786 {
02787 int32_t value;
02788 while (((value=fgetc(file)) != EOF) && (value!=']'))
02789 {
02790 if (value=='\n')
02791 line++;
02792 }
02793
02794 if (value==EOF)
02795 error(line, "expected \"]\" in input file");
02796 }
02797
02798 bool CHMM::comma_or_space(FILE* file)
02799 {
02800 int32_t value;
02801 while (((value=fgetc(file)) != EOF) && (value!=',') && (value!=';') && (value!=']'))
02802 {
02803 if (value=='\n')
02804 line++;
02805 }
02806 if (value==']')
02807 {
02808 ungetc(value, file);
02809 SG_ERROR( "found ']' instead of ';' or ','\n") ;
02810 return false ;
02811 } ;
02812
02813 if (value==EOF)
02814 error(line, "expected \";\" or \",\" in input file");
02815
02816 while (((value=fgetc(file)) != EOF) && (isspace(value)))
02817 {
02818 if (value=='\n')
02819 line++;
02820 }
02821 ungetc(value, file);
02822 return true ;
02823 }
02824
02825 bool CHMM::get_numbuffer(FILE* file, char* buffer, int32_t length)
02826 {
02827 signed char value;
02828
02829 while (((value=fgetc(file)) != EOF) &&
02830 !isdigit(value) && (value!='A')
02831 && (value!='C') && (value!='G') && (value!='T')
02832 && (value!='N') && (value!='n')
02833 && (value!='.') && (value!='-') && (value!='e') && (value!=']'))
02834 {
02835 if (value=='\n')
02836 line++;
02837 }
02838 if (value==']')
02839 {
02840 ungetc(value,file) ;
02841 return false ;
02842 } ;
02843 if (value!=EOF)
02844 {
02845 int32_t i=0;
02846 switch (value)
02847 {
02848 case 'A':
02849 value='0' +CAlphabet::B_A;
02850 break;
02851 case 'C':
02852 value='0' +CAlphabet::B_C;
02853 break;
02854 case 'G':
02855 value='0' +CAlphabet::B_G;
02856 break;
02857 case 'T':
02858 value='0' +CAlphabet::B_T;
02859 break;
02860 };
02861
02862 buffer[i++]=value;
02863
02864 while (((value=fgetc(file)) != EOF) &&
02865 (isdigit(value) || (value=='.') || (value=='-') || (value=='e')
02866 || (value=='A') || (value=='C') || (value=='G')|| (value=='T')
02867 || (value=='N') || (value=='n')) && (i<length))
02868 {
02869 switch (value)
02870 {
02871 case 'A':
02872 value='0' +CAlphabet::B_A;
02873 break;
02874 case 'C':
02875 value='0' +CAlphabet::B_C;
02876 break;
02877 case 'G':
02878 value='0' +CAlphabet::B_G;
02879 break;
02880 case 'T':
02881 value='0' +CAlphabet::B_T;
02882 break;
02883 case '1': case '2': case'3': case '4': case'5':
02884 case '6': case '7': case'8': case '9': case '0': break ;
02885 case '.': case 'e': case '-': break ;
02886 default:
02887 SG_ERROR( "found crap: %i %c (pos:%li)\n",i,value,ftell(file)) ;
02888 };
02889 buffer[i++]=value;
02890 }
02891 ungetc(value, file);
02892 buffer[i]='\0';
02893
02894 return (i<=length) && (i>0);
02895 }
02896 return false;
02897 }
02898
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909
02910
02911
02912
02913
02914
02915
02916
02917
02918
02919
02920
02921
02922
02923
02924
02925
02926
02927
02928
02929
02930
02931
02932
02933
02934 bool CHMM::load_model(FILE* file)
02935 {
02936 int32_t received_params=0;
02937
02938 bool result=false;
02939 E_STATE state=INITIAL;
02940 char buffer[1024];
02941
02942 line=1;
02943 int32_t i,j;
02944
02945 if (file)
02946 {
02947 while (state!=END)
02948 {
02949 int32_t value=fgetc(file);
02950
02951 if (value=='\n')
02952 line++;
02953 if (value==EOF)
02954 state=END;
02955
02956 switch (state)
02957 {
02958 case INITIAL:
02959 if (value=='N')
02960 {
02961 if (received_params & GOTN)
02962 error(line, "in model file: \"p double defined\"");
02963 else
02964 state=GET_N;
02965 }
02966 else if (value=='M')
02967 {
02968 if (received_params & GOTM)
02969 error(line, "in model file: \"p double defined\"");
02970 else
02971 state=GET_M;
02972 }
02973 else if (value=='%')
02974 {
02975 state=COMMENT;
02976 }
02977 case ARRAYs:
02978 if (value=='p')
02979 {
02980 if (received_params & GOTp)
02981 error(line, "in model file: \"p double defined\"");
02982 else
02983 state=GET_p;
02984 }
02985 if (value=='q')
02986 {
02987 if (received_params & GOTq)
02988 error(line, "in model file: \"q double defined\"");
02989 else
02990 state=GET_q;
02991 }
02992 else if (value=='a')
02993 {
02994 if (received_params & GOTa)
02995 error(line, "in model file: \"a double defined\"");
02996 else
02997 state=GET_a;
02998 }
02999 else if (value=='b')
03000 {
03001 if (received_params & GOTb)
03002 error(line, "in model file: \"b double defined\"");
03003 else
03004 state=GET_b;
03005 }
03006 else if (value=='%')
03007 {
03008 state=COMMENT;
03009 }
03010 break;
03011 case GET_N:
03012 if (value=='=')
03013 {
03014 if (get_numbuffer(file, buffer, 4))
03015 {
03016 this->N= atoi(buffer);
03017 received_params|=GOTN;
03018 state= (received_params == (GOTN | GOTM | GOTO)) ? ARRAYs : INITIAL;
03019 }
03020 else
03021 state=END;
03022 }
03023 break;
03024 case GET_M:
03025 if (value=='=')
03026 {
03027 if (get_numbuffer(file, buffer, 4))
03028 {
03029 this->M= atoi(buffer);
03030 received_params|=GOTM;
03031 state= (received_params == (GOTN | GOTM | GOTO)) ? ARRAYs : INITIAL;
03032 }
03033 else
03034 state=END;
03035 }
03036 break;
03037 case GET_a:
03038 if (value=='=')
03039 {
03040 float64_t f;
03041
03042 transition_matrix_a=SG_MALLOC(float64_t, N*N);
03043 open_bracket(file);
03044 for (i=0; i<this->N; i++)
03045 {
03046 open_bracket(file);
03047
03048 for (j=0; j<this->N ; j++)
03049 {
03050
03051 if (fscanf( file, "%le", &f ) != 1)
03052 error(line, "float64_t expected");
03053 else
03054 set_a(i,j, f);
03055
03056 if (j<this->N-1)
03057 comma_or_space(file);
03058 else
03059 close_bracket(file);
03060 }
03061
03062 if (i<this->N-1)
03063 comma_or_space(file);
03064 else
03065 close_bracket(file);
03066 }
03067 received_params|=GOTa;
03068 }
03069 state= (received_params == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
03070 break;
03071 case GET_b:
03072 if (value=='=')
03073 {
03074 float64_t f;
03075
03076 observation_matrix_b=SG_MALLOC(float64_t, N*M);
03077 open_bracket(file);
03078 for (i=0; i<this->N; i++)
03079 {
03080 open_bracket(file);
03081
03082 for (j=0; j<this->M ; j++)
03083 {
03084
03085 if (fscanf( file, "%le", &f ) != 1)
03086 error(line, "float64_t expected");
03087 else
03088 set_b(i,j, f);
03089
03090 if (j<this->M-1)
03091 comma_or_space(file);
03092 else
03093 close_bracket(file);
03094 }
03095
03096 if (i<this->N-1)
03097 comma_or_space(file);
03098 else
03099 close_bracket(file);
03100 }
03101 received_params|=GOTb;
03102 }
03103 state= ((received_params & (GOTa | GOTb | GOTp | GOTq)) == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
03104 break;
03105 case GET_p:
03106 if (value=='=')
03107 {
03108 float64_t f;
03109
03110 initial_state_distribution_p=SG_MALLOC(float64_t, N);
03111 open_bracket(file);
03112 for (i=0; i<this->N ; i++)
03113 {
03114 if (fscanf( file, "%le", &f ) != 1)
03115 error(line, "float64_t expected");
03116 else
03117 set_p(i, f);
03118
03119 if (i<this->N-1)
03120 comma_or_space(file);
03121 else
03122 close_bracket(file);
03123 }
03124 received_params|=GOTp;
03125 }
03126 state= (received_params == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
03127 break;
03128 case GET_q:
03129 if (value=='=')
03130 {
03131 float64_t f;
03132
03133 end_state_distribution_q=SG_MALLOC(float64_t, N);
03134 open_bracket(file);
03135 for (i=0; i<this->N ; i++)
03136 {
03137 if (fscanf( file, "%le", &f ) != 1)
03138 error(line, "float64_t expected");
03139 else
03140 set_q(i, f);
03141
03142 if (i<this->N-1)
03143 comma_or_space(file);
03144 else
03145 close_bracket(file);
03146 }
03147 received_params|=GOTq;
03148 }
03149 state= (received_params == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
03150 break;
03151 case COMMENT:
03152 if (value==EOF)
03153 state=END;
03154 else if (value=='\n')
03155 {
03156 line++;
03157 state=INITIAL;
03158 }
03159 break;
03160
03161 default:
03162 break;
03163 }
03164 }
03165 result= (received_params== (GOTa | GOTb | GOTp | GOTq | GOTN | GOTM | GOTO));
03166 }
03167
03168 SG_WARNING( "not normalizing anymore, call normalize_hmm to make sure the hmm is valid!!\n");
03170 return result;
03171 }
03172
03173
03174
03175
03176
03177
03178
03179
03180
03181
03182
03183
03184
03185
03186
03187
03188
03189
03190
03191
03192
03193
03194
03195
03196
03197
03198
03199
03200
03201
03202
03203
03204
03205
03206
03207
03208
03209
03210
03211
03212
03213
03214
03215
03216
03217
03218
03219
03220
03221
03222
03223
03224
03225
03226
03227
03228
03229
03230
03231
03232 bool CHMM::load_definitions(FILE* file, bool verbose, bool _initialize)
03233 {
03234 if (model)
03235 delete model ;
03236 model=new Model();
03237
03238 int32_t received_params=0x0000000;
03239 char buffer[1024];
03240
03241 bool result=false;
03242 E_STATE state=INITIAL;
03243
03244 {
03245 model->set_learn_a(0, -1);
03246 model->set_learn_a(1, -1);
03247 model->set_const_a(0, -1);
03248 model->set_const_a(1, -1);
03249 model->set_const_a_val(0, 1.0);
03250 model->set_learn_b(0, -1);
03251 model->set_const_b(0, -1);
03252 model->set_const_b_val(0, 1.0);
03253 model->set_learn_p(0, -1);
03254 model->set_learn_q(0, -1);
03255 model->set_const_p(0, -1);
03256 model->set_const_q(0, -1);
03257 } ;
03258
03259 line=1;
03260
03261 if (file)
03262 {
03263 while (state!=END)
03264 {
03265 int32_t value=fgetc(file);
03266
03267 if (value=='\n')
03268 line++;
03269
03270 if (value==EOF)
03271 state=END;
03272
03273 switch (state)
03274 {
03275 case INITIAL:
03276 if (value=='l')
03277 {
03278 if (fgetc(file)=='e' && fgetc(file)=='a' && fgetc(file)=='r' && fgetc(file)=='n' && fgetc(file)=='_')
03279 {
03280 switch(fgetc(file))
03281 {
03282 case 'a':
03283 state=GET_learn_a;
03284 break;
03285 case 'b':
03286 state=GET_learn_b;
03287 break;
03288 case 'p':
03289 state=GET_learn_p;
03290 break;
03291 case 'q':
03292 state=GET_learn_q;
03293 break;
03294 default:
03295 error(line, "a,b,p or q expected in train definition file");
03296 };
03297 }
03298 }
03299 else if (value=='c')
03300 {
03301 if (fgetc(file)=='o' && fgetc(file)=='n' && fgetc(file)=='s'
03302 && fgetc(file)=='t' && fgetc(file)=='_')
03303 {
03304 switch(fgetc(file))
03305 {
03306 case 'a':
03307 state=GET_const_a;
03308 break;
03309 case 'b':
03310 state=GET_const_b;
03311 break;
03312 case 'p':
03313 state=GET_const_p;
03314 break;
03315 case 'q':
03316 state=GET_const_q;
03317 break;
03318 default:
03319 error(line, "a,b,p or q expected in train definition file");
03320 };
03321 }
03322 }
03323 else if (value=='%')
03324 {
03325 state=COMMENT;
03326 }
03327 else if (value==EOF)
03328 {
03329 state=END;
03330 }
03331 break;
03332 case GET_learn_a:
03333 if (value=='=')
03334 {
03335 open_bracket(file);
03336 bool finished=false;
03337 int32_t i=0;
03338
03339 if (verbose)
03340 SG_DEBUG( "\nlearn for transition matrix: ") ;
03341 while (!finished)
03342 {
03343 open_bracket(file);
03344
03345 if (get_numbuffer(file, buffer, 4))
03346 {
03347 value=atoi(buffer);
03348 model->set_learn_a(i++, value);
03349
03350 if (value<0)
03351 {
03352 finished=true;
03353 break;
03354 }
03355 if (value>=N)
03356 SG_ERROR( "invalid value for learn_a(%i,0): %i\n",i/2,(int)value) ;
03357 }
03358 else
03359 break;
03360
03361 comma_or_space(file);
03362
03363 if (get_numbuffer(file, buffer, 4))
03364 {
03365 value=atoi(buffer);
03366 model->set_learn_a(i++, value);
03367
03368 if (value<0)
03369 {
03370 finished=true;
03371 break;
03372 }
03373 if (value>=N)
03374 SG_ERROR( "invalid value for learn_a(%i,1): %i\n",i/2-1,(int)value) ;
03375
03376 }
03377 else
03378 break;
03379 close_bracket(file);
03380 }
03381 close_bracket(file);
03382 if (verbose)
03383 SG_DEBUG( "%i Entries",(int)(i/2)) ;
03384
03385 if (finished)
03386 {
03387 received_params|=GOTlearn_a;
03388
03389 state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
03390 }
03391 else
03392 state=END;
03393 }
03394 break;
03395 case GET_learn_b:
03396 if (value=='=')
03397 {
03398 open_bracket(file);
03399 bool finished=false;
03400 int32_t i=0;
03401
03402 if (verbose)
03403 SG_DEBUG( "\nlearn for emission matrix: ") ;
03404
03405 while (!finished)
03406 {
03407 open_bracket(file);
03408
03409 int32_t combine=0;
03410
03411 for (int32_t j=0; j<2; j++)
03412 {
03413 if (get_numbuffer(file, buffer, 4))
03414 {
03415 value=atoi(buffer);
03416
03417 if (j==0)
03418 {
03419 model->set_learn_b(i++, value);
03420
03421 if (value<0)
03422 {
03423 finished=true;
03424 break;
03425 }
03426 if (value>=N)
03427 SG_ERROR( "invalid value for learn_b(%i,0): %i\n",i/2,(int)value) ;
03428 }
03429 else
03430 combine=value;
03431 }
03432 else
03433 break;
03434
03435 if (j<1)
03436 comma_or_space(file);
03437 else
03438 close_bracket(file);
03439 }
03440 model->set_learn_b(i++, combine);
03441 if (combine>=M)
03442
03443 SG_ERROR("invalid value for learn_b(%i,1): %i\n",i/2-1,(int)value) ;
03444 }
03445 close_bracket(file);
03446 if (verbose)
03447 SG_DEBUG( "%i Entries",(int)(i/2-1)) ;
03448
03449 if (finished)
03450 {
03451 received_params|=GOTlearn_b;
03452 state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
03453 }
03454 else
03455 state=END;
03456 }
03457 break;
03458 case GET_learn_p:
03459 if (value=='=')
03460 {
03461 open_bracket(file);
03462 bool finished=false;
03463 int32_t i=0;
03464
03465 if (verbose)
03466 SG_DEBUG( "\nlearn start states: ") ;
03467 while (!finished)
03468 {
03469 if (get_numbuffer(file, buffer, 4))
03470 {
03471 value=atoi(buffer);
03472
03473 model->set_learn_p(i++, value);
03474
03475 if (value<0)
03476 {
03477 finished=true;
03478 break;
03479 }
03480 if (value>=N)
03481 SG_ERROR( "invalid value for learn_p(%i): %i\n",i-1,(int)value) ;
03482 }
03483 else
03484 break;
03485
03486 comma_or_space(file);
03487 }
03488
03489 close_bracket(file);
03490 if (verbose)
03491 SG_DEBUG( "%i Entries",i-1) ;
03492
03493 if (finished)
03494 {
03495 received_params|=GOTlearn_p;
03496 state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
03497 }
03498 else
03499 state=END;
03500 }
03501 break;
03502 case GET_learn_q:
03503 if (value=='=')
03504 {
03505 open_bracket(file);
03506 bool finished=false;
03507 int32_t i=0;
03508
03509 if (verbose)
03510 SG_DEBUG( "\nlearn terminal states: ") ;
03511 while (!finished)
03512 {
03513 if (get_numbuffer(file, buffer, 4))
03514 {
03515 value=atoi(buffer);
03516 model->set_learn_q(i++, value);
03517
03518 if (value<0)
03519 {
03520 finished=true;
03521 break;
03522 }
03523 if (value>=N)
03524 SG_ERROR( "invalid value for learn_q(%i): %i\n",i-1,(int)value) ;
03525 }
03526 else
03527 break;
03528
03529 comma_or_space(file);
03530 }
03531
03532 close_bracket(file);
03533 if (verbose)
03534 SG_DEBUG( "%i Entries",i-1) ;
03535
03536 if (finished)
03537 {
03538 received_params|=GOTlearn_q;
03539 state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
03540 }
03541 else
03542 state=END;
03543 }
03544 break;
03545 case GET_const_a:
03546 if (value=='=')
03547 {
03548 open_bracket(file);
03549 bool finished=false;
03550 int32_t i=0;
03551
03552 if (verbose)
03553 #ifdef USE_HMMDEBUG
03554 SG_DEBUG( "\nconst for transition matrix: \n") ;
03555 #else
03556 SG_DEBUG( "\nconst for transition matrix: ") ;
03557 #endif
03558 while (!finished)
03559 {
03560 open_bracket(file);
03561
03562 if (get_numbuffer(file, buffer, 4))
03563 {
03564 value=atoi(buffer);
03565 model->set_const_a(i++, value);
03566
03567 if (value<0)
03568 {
03569 finished=true;
03570 model->set_const_a(i++, value);
03571 model->set_const_a_val((int32_t)i/2 - 1, value);
03572 break;
03573 }
03574 if (value>=N)
03575 SG_ERROR( "invalid value for const_a(%i,0): %i\n",i/2,(int)value) ;
03576 }
03577 else
03578 break;
03579
03580 comma_or_space(file);
03581
03582 if (get_numbuffer(file, buffer, 4))
03583 {
03584 value=atoi(buffer);
03585 model->set_const_a(i++, value);
03586
03587 if (value<0)
03588 {
03589 finished=true;
03590 model->set_const_a_val((int32_t)i/2 - 1, value);
03591 break;
03592 }
03593 if (value>=N)
03594 SG_ERROR( "invalid value for const_a(%i,1): %i\n",i/2-1,(int)value) ;
03595 }
03596 else
03597 break;
03598
03599 if (!comma_or_space(file))
03600 model->set_const_a_val((int32_t)i/2 - 1, 1.0);
03601 else
03602 if (get_numbuffer(file, buffer, 10))
03603 {
03604 float64_t dvalue=atof(buffer);
03605 model->set_const_a_val((int32_t)i/2 - 1, dvalue);
03606 if (dvalue<0)
03607 {
03608 finished=true;
03609 break;
03610 }
03611 if ((dvalue>1.0) || (dvalue<0.0))
03612 SG_ERROR( "invalid value for const_a_val(%i): %e\n",(int)i/2-1,dvalue) ;
03613 }
03614 else
03615 model->set_const_a_val((int32_t)i/2 - 1, 1.0);
03616
03617 #ifdef USE_HMMDEBUG
03618 if (verbose)
03619 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)) ;
03620 #endif
03621 close_bracket(file);
03622 }
03623 close_bracket(file);
03624 if (verbose)
03625 SG_DEBUG( "%i Entries",(int)i/2-1) ;
03626
03627 if (finished)
03628 {
03629 received_params|=GOTconst_a;
03630 state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
03631 }
03632 else
03633 state=END;
03634 }
03635 break;
03636
03637 case GET_const_b:
03638 if (value=='=')
03639 {
03640 open_bracket(file);
03641 bool finished=false;
03642 int32_t i=0;
03643
03644 if (verbose)
03645 #ifdef USE_HMMDEBUG
03646 SG_DEBUG( "\nconst for emission matrix: \n") ;
03647 #else
03648 SG_DEBUG( "\nconst for emission matrix: ") ;
03649 #endif
03650 while (!finished)
03651 {
03652 open_bracket(file);
03653 int32_t combine=0;
03654 for (int32_t j=0; j<3; j++)
03655 {
03656 if (get_numbuffer(file, buffer, 10))
03657 {
03658 if (j==0)
03659 {
03660 value=atoi(buffer);
03661
03662 model->set_const_b(i++, value);
03663
03664 if (value<0)
03665 {
03666 finished=true;
03667
03668 break;
03669 }
03670 if (value>=N)
03671 SG_ERROR( "invalid value for const_b(%i,0): %i\n",i/2-1,(int)value) ;
03672 }
03673 else if (j==2)
03674 {
03675 float64_t dvalue=atof(buffer);
03676 model->set_const_b_val((int32_t)(i-1)/2, dvalue);
03677 if (dvalue<0)
03678 {
03679 finished=true;
03680 break;
03681 } ;
03682 if ((dvalue>1.0) || (dvalue<0.0))
03683 SG_ERROR( "invalid value for const_b_val(%i,1): %e\n",i/2-1,dvalue) ;
03684 }
03685 else
03686 {
03687 value=atoi(buffer);
03688 combine= value;
03689 } ;
03690 }
03691 else
03692 {
03693 if (j==2)
03694 model->set_const_b_val((int32_t)(i-1)/2, 1.0);
03695 break;
03696 } ;
03697 if (j<2)
03698 if ((!comma_or_space(file)) && (j==1))
03699 {
03700 model->set_const_b_val((int32_t)(i-1)/2, 1.0) ;
03701 break ;
03702 } ;
03703 }
03704 close_bracket(file);
03705 model->set_const_b(i++, combine);
03706 if (combine>=M)
03707 SG_ERROR("invalid value for const_b(%i,1): %i\n",i/2-1, combine) ;
03708 #ifdef USE_HMMDEBUG
03709 if (verbose && !finished)
03710 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)) ;
03711 #endif
03712 }
03713 close_bracket(file);
03714 if (verbose)
03715 SG_ERROR( "%i Entries",(int)i/2-1) ;
03716
03717 if (finished)
03718 {
03719 received_params|=GOTconst_b;
03720 state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
03721 }
03722 else
03723 state=END;
03724 }
03725 break;
03726 case GET_const_p:
03727 if (value=='=')
03728 {
03729 open_bracket(file);
03730 bool finished=false;
03731 int32_t i=0;
03732
03733 if (verbose)
03734 #ifdef USE_HMMDEBUG
03735 SG_DEBUG( "\nconst for start states: \n") ;
03736 #else
03737 SG_DEBUG( "\nconst for start states: ") ;
03738 #endif
03739 while (!finished)
03740 {
03741 open_bracket(file);
03742
03743 if (get_numbuffer(file, buffer, 4))
03744 {
03745 value=atoi(buffer);
03746 model->set_const_p(i, value);
03747
03748 if (value<0)
03749 {
03750 finished=true;
03751 model->set_const_p_val(i++, value);
03752 break;
03753 }
03754 if (value>=N)
03755 SG_ERROR( "invalid value for const_p(%i): %i\n",i,(int)value) ;
03756
03757 }
03758 else
03759 break;
03760
03761 if (!comma_or_space(file))
03762 model->set_const_p_val(i++, 1.0);
03763 else
03764 if (get_numbuffer(file, buffer, 10))
03765 {
03766 float64_t dvalue=atof(buffer);
03767 model->set_const_p_val(i++, dvalue);
03768 if (dvalue<0)
03769 {
03770 finished=true;
03771 break;
03772 }
03773 if ((dvalue>1) || (dvalue<0))
03774 SG_ERROR( "invalid value for const_p_val(%i): %e\n",i,dvalue) ;
03775 }
03776 else
03777 model->set_const_p_val(i++, 1.0);
03778
03779 close_bracket(file);
03780
03781 #ifdef USE_HMMDEBUG
03782 if (verbose)
03783 SG_DEBUG("const_p(%i)=%e\n", model->get_const_p(i-1),model->get_const_p_val(i-1)) ;
03784 #endif
03785 }
03786 if (verbose)
03787 SG_DEBUG( "%i Entries",i-1) ;
03788
03789 close_bracket(file);
03790
03791 if (finished)
03792 {
03793 received_params|=GOTconst_p;
03794 state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
03795 }
03796 else
03797 state=END;
03798 }
03799 break;
03800 case GET_const_q:
03801 if (value=='=')
03802 {
03803 open_bracket(file);
03804 bool finished=false;
03805 if (verbose)
03806 #ifdef USE_HMMDEBUG
03807 SG_DEBUG( "\nconst for terminal states: \n") ;
03808 #else
03809 SG_DEBUG( "\nconst for terminal states: ") ;
03810 #endif
03811 int32_t i=0;
03812
03813 while (!finished)
03814 {
03815 open_bracket(file) ;
03816 if (get_numbuffer(file, buffer, 4))
03817 {
03818 value=atoi(buffer);
03819 model->set_const_q(i, value);
03820 if (value<0)
03821 {
03822 finished=true;
03823 model->set_const_q_val(i++, value);
03824 break;
03825 }
03826 if (value>=N)
03827 SG_ERROR( "invalid value for const_q(%i): %i\n",i,(int)value) ;
03828 }
03829 else
03830 break;
03831
03832 if (!comma_or_space(file))
03833 model->set_const_q_val(i++, 1.0);
03834 else
03835 if (get_numbuffer(file, buffer, 10))
03836 {
03837 float64_t dvalue=atof(buffer);
03838 model->set_const_q_val(i++, dvalue);
03839 if (dvalue<0)
03840 {
03841 finished=true;
03842 break;
03843 }
03844 if ((dvalue>1) || (dvalue<0))
03845 SG_ERROR("invalid value for const_q_val(%i): %e\n",i,(double) dvalue) ;
03846 }
03847 else
03848 model->set_const_q_val(i++, 1.0);
03849
03850 close_bracket(file);
03851 #ifdef USE_HMMDEBUG
03852 if (verbose)
03853 SG_DEBUG("const_q(%i)=%e\n", model->get_const_q(i-1),model->get_const_q_val(i-1)) ;
03854 #endif
03855 }
03856 if (verbose)
03857 SG_DEBUG( "%i Entries",i-1) ;
03858
03859 close_bracket(file);
03860
03861 if (finished)
03862 {
03863 received_params|=GOTconst_q;
03864 state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
03865 }
03866 else
03867 state=END;
03868 }
03869 break;
03870 case COMMENT:
03871 if (value==EOF)
03872 state=END;
03873 else if (value=='\n')
03874 state=INITIAL;
03875 break;
03876
03877 default:
03878 break;
03879 }
03880 }
03881 }
03882
03883
03884
03885
03886
03887 result=1 ;
03888 if (result)
03889 {
03890 model->sort_learn_a() ;
03891 model->sort_learn_b() ;
03892 if (_initialize)
03893 {
03894 init_model_defined(); ;
03895 convert_to_log();
03896 } ;
03897 }
03898 if (verbose)
03899 SG_DEBUG( "\n") ;
03900 return result;
03901 }
03902
03903
03904
03905
03906
03907
03908
03909
03910
03911
03912
03913
03914
03915
03916
03917
03918
03919
03920
03921
03922
03923
03924
03925
03926
03927
03928
03929
03930
03931
03932
03933
03934
03935
03936
03937 bool CHMM::save_model(FILE* file)
03938 {
03939 bool result=false;
03940 int32_t i,j;
03941 const float32_t NAN_REPLACEMENT = (float32_t) CMath::ALMOST_NEG_INFTY ;
03942
03943 if (file)
03944 {
03945 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");
03946 fprintf(file,"N=%d;\n",N);
03947 fprintf(file,"M=%d;\n",M);
03948
03949 fprintf(file,"p=[");
03950 for (i=0; i<N; i++)
03951 {
03952 if (i<N-1) {
03953 if (CMath::is_finite(get_p(i)))
03954 fprintf(file, "%e,", (double)get_p(i));
03955 else
03956 fprintf(file, "%f,", NAN_REPLACEMENT);
03957 }
03958 else {
03959 if (CMath::is_finite(get_p(i)))
03960 fprintf(file, "%e", (double)get_p(i));
03961 else
03962 fprintf(file, "%f", NAN_REPLACEMENT);
03963 }
03964 }
03965
03966 fprintf(file,"];\n\nq=[");
03967 for (i=0; i<N; i++)
03968 {
03969 if (i<N-1) {
03970 if (CMath::is_finite(get_q(i)))
03971 fprintf(file, "%e,", (double)get_q(i));
03972 else
03973 fprintf(file, "%f,", NAN_REPLACEMENT);
03974 }
03975 else {
03976 if (CMath::is_finite(get_q(i)))
03977 fprintf(file, "%e", (double)get_q(i));
03978 else
03979 fprintf(file, "%f", NAN_REPLACEMENT);
03980 }
03981 }
03982 fprintf(file,"];\n\na=[");
03983
03984 for (i=0; i<N; i++)
03985 {
03986 fprintf(file, "\t[");
03987
03988 for (j=0; j<N; j++)
03989 {
03990 if (j<N-1) {
03991 if (CMath::is_finite(get_a(i,j)))
03992 fprintf(file, "%e,", (double)get_a(i,j));
03993 else
03994 fprintf(file, "%f,", NAN_REPLACEMENT);
03995 }
03996 else {
03997 if (CMath::is_finite(get_a(i,j)))
03998 fprintf(file, "%e];\n", (double)get_a(i,j));
03999 else
04000 fprintf(file, "%f];\n", NAN_REPLACEMENT);
04001 }
04002 }
04003 }
04004
04005 fprintf(file," ];\n\nb=[");
04006
04007 for (i=0; i<N; i++)
04008 {
04009 fprintf(file, "\t[");
04010
04011 for (j=0; j<M; j++)
04012 {
04013 if (j<M-1) {
04014 if (CMath::is_finite(get_b(i,j)))
04015 fprintf(file, "%e,", (double)get_b(i,j));
04016 else
04017 fprintf(file, "%f,", NAN_REPLACEMENT);
04018 }
04019 else {
04020 if (CMath::is_finite(get_b(i,j)))
04021 fprintf(file, "%e];\n", (double)get_b(i,j));
04022 else
04023 fprintf(file, "%f];\n", NAN_REPLACEMENT);
04024 }
04025 }
04026 }
04027 result= (fprintf(file," ];\n") == 5);
04028 }
04029
04030 return result;
04031 }
04032
04033 T_STATES* CHMM::get_path(int32_t dim, float64_t& prob)
04034 {
04035 T_STATES* result = NULL;
04036
04037 prob = best_path(dim);
04038 result = SG_MALLOC(T_STATES, p_observations->get_vector_length(dim));
04039
04040 for (int32_t i=0; i<p_observations->get_vector_length(dim); i++)
04041 result[i]=PATH(dim)[i];
04042
04043 return result;
04044 }
04045
04046 bool CHMM::save_path(FILE* file)
04047 {
04048 bool result=false;
04049
04050 if (file)
04051 {
04052 for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
04053 {
04054 if (dim%100==0)
04055 SG_PRINT( "%i..", dim) ;
04056 float64_t prob = best_path(dim);
04057 fprintf(file,"%i. path probability:%e\nstate sequence:\n", dim, prob);
04058 for (int32_t i=0; i<p_observations->get_vector_length(dim)-1; i++)
04059 fprintf(file,"%d ", PATH(dim)[i]);
04060 fprintf(file,"%d", PATH(dim)[p_observations->get_vector_length(dim)-1]);
04061 fprintf(file,"\n\n") ;
04062 }
04063 SG_DONE();
04064 result=true;
04065 }
04066
04067 return result;
04068 }
04069
04070 bool CHMM::save_likelihood_bin(FILE* file)
04071 {
04072 bool result=false;
04073
04074 if (file)
04075 {
04076 for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
04077 {
04078 float32_t prob= (float32_t) model_probability(dim);
04079 fwrite(&prob, sizeof(float32_t), 1, file);
04080 }
04081 result=true;
04082 }
04083
04084 return result;
04085 }
04086
04087 bool CHMM::save_likelihood(FILE* file)
04088 {
04089 bool result=false;
04090
04091 if (file)
04092 {
04093 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");
04094
04095 fprintf(file, "P=[");
04096 for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
04097 fprintf(file, "%e ", (double) model_probability(dim));
04098
04099 fprintf(file,"];");
04100 result=true;
04101 }
04102
04103 return result;
04104 }
04105
04106 #define FLOATWRITE(file, value) { float32_t rrr=float32_t(value); fwrite(&rrr, sizeof(float32_t), 1, file); num_floats++;}
04107
04108 bool CHMM::save_model_bin(FILE* file)
04109 {
04110 int32_t i,j,q, num_floats=0 ;
04111 if (!model)
04112 {
04113 if (file)
04114 {
04115
04116 FLOATWRITE(file, (float32_t)CMath::INFTY);
04117 FLOATWRITE(file, (float32_t) 1);
04118
04119
04120 for (i=0; i<N; i++)
04121 FLOATWRITE(file, get_p(i));
04122 SG_INFO( "wrote %i parameters for p\n",N) ;
04123
04124 for (i=0; i<N; i++)
04125 FLOATWRITE(file, get_q(i)) ;
04126 SG_INFO( "wrote %i parameters for q\n",N) ;
04127
04128
04129 for (i=0; i<N; i++)
04130 for (j=0; j<N; j++)
04131 FLOATWRITE(file, get_a(i,j));
04132 SG_INFO( "wrote %i parameters for a\n",N*N) ;
04133
04134 for (i=0; i<N; i++)
04135 for (j=0; j<M; j++)
04136 FLOATWRITE(file, get_b(i,j));
04137 SG_INFO( "wrote %i parameters for b\n",N*M) ;
04138
04139
04140 FLOATWRITE(file, (float32_t)CMath::INFTY);
04141 FLOATWRITE(file, (float32_t) 3);
04142
04143
04144 FLOATWRITE(file, (float32_t) N);
04145 FLOATWRITE(file, (float32_t) N);
04146 FLOATWRITE(file, (float32_t) N*N);
04147 FLOATWRITE(file, (float32_t) N*M);
04148 FLOATWRITE(file, (float32_t) N);
04149 FLOATWRITE(file, (float32_t) M);
04150 } ;
04151 }
04152 else
04153 {
04154 if (file)
04155 {
04156 int32_t num_p, num_q, num_a, num_b ;
04157
04158 FLOATWRITE(file, (float32_t)CMath::INFTY);
04159 FLOATWRITE(file, (float32_t) 2);
04160
04161 for (i=0; model->get_learn_p(i)>=0; i++)
04162 FLOATWRITE(file, get_p(model->get_learn_p(i)));
04163 num_p=i ;
04164 SG_INFO( "wrote %i parameters for p\n",num_p) ;
04165
04166 for (i=0; model->get_learn_q(i)>=0; i++)
04167 FLOATWRITE(file, get_q(model->get_learn_q(i)));
04168 num_q=i ;
04169 SG_INFO( "wrote %i parameters for q\n",num_q) ;
04170
04171
04172 for (q=0; model->get_learn_a(q,1)>=0; q++)
04173 {
04174 i=model->get_learn_a(q,0) ;
04175 j=model->get_learn_a(q,1) ;
04176 FLOATWRITE(file, (float32_t)i);
04177 FLOATWRITE(file, (float32_t)j);
04178 FLOATWRITE(file, get_a(i,j));
04179 } ;
04180 num_a=q ;
04181 SG_INFO( "wrote %i parameters for a\n",num_a) ;
04182
04183 for (q=0; model->get_learn_b(q,0)>=0; q++)
04184 {
04185 i=model->get_learn_b(q,0) ;
04186 j=model->get_learn_b(q,1) ;
04187 FLOATWRITE(file, (float32_t)i);
04188 FLOATWRITE(file, (float32_t)j);
04189 FLOATWRITE(file, get_b(i,j));
04190 } ;
04191 num_b=q ;
04192 SG_INFO( "wrote %i parameters for b\n",num_b) ;
04193
04194
04195 FLOATWRITE(file, (float32_t)CMath::INFTY);
04196 FLOATWRITE(file, (float32_t) 3);
04197
04198
04199 FLOATWRITE(file, (float32_t) num_p);
04200 FLOATWRITE(file, (float32_t) num_q);
04201 FLOATWRITE(file, (float32_t) num_a);
04202 FLOATWRITE(file, (float32_t) num_b);
04203 FLOATWRITE(file, (float32_t) N);
04204 FLOATWRITE(file, (float32_t) M);
04205 } ;
04206 } ;
04207 return true ;
04208 }
04209
04210 bool CHMM::save_path_derivatives(FILE* logfile)
04211 {
04212 int32_t dim,i,j;
04213
04214 if (logfile)
04215 {
04216 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");
04217 fprintf(logfile,"%% dPr[...]=[ [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
04218 fprintf(logfile,"%% [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
04219 fprintf(logfile,"%% ............................. \n");
04220 fprintf(logfile,"%% [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_MM]\n");
04221 fprintf(logfile,"%% ];\n%%\n\ndPr(log()) = [\n");
04222 }
04223 else
04224 return false ;
04225
04226 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
04227 {
04228 best_path(dim);
04229
04230 fprintf(logfile, "[ ");
04231
04232
04233 for (i=0; i<N; i++)
04234 fprintf(logfile,"%e, ", (double) path_derivative_p(i,dim) );
04235
04236 for (i=0; i<N; i++)
04237 fprintf(logfile,"%e, ", (double) path_derivative_q(i,dim) );
04238
04239
04240 for (i=0; i<N; i++)
04241 for (j=0; j<N; j++)
04242 fprintf(logfile, "%e,", (double) path_derivative_a(i,j,dim) );
04243
04244 for (i=0; i<N; i++)
04245 for (j=0; j<M; j++)
04246 fprintf(logfile, "%e,", (double) path_derivative_b(i,j,dim) );
04247
04248 fseek(logfile,ftell(logfile)-1,SEEK_SET);
04249 fprintf(logfile, " ];\n");
04250 }
04251
04252 fprintf(logfile, "];");
04253
04254 return true ;
04255 }
04256
04257 bool CHMM::save_path_derivatives_bin(FILE* logfile)
04258 {
04259 bool result=false;
04260 int32_t dim,i,j,q;
04261 float64_t prob=0 ;
04262 int32_t num_floats=0 ;
04263
04264 float64_t sum_prob=0.0 ;
04265 if (!model)
04266 SG_WARNING( "No definitions loaded -- writing derivatives of all weights\n") ;
04267 else
04268 SG_INFO( "writing derivatives of changed weights only\n") ;
04269
04270 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
04271 {
04272 if (dim%100==0)
04273 {
04274 SG_PRINT( ".") ;
04275
04276 } ;
04277
04278 prob=best_path(dim);
04279 sum_prob+=prob ;
04280
04281 if (!model)
04282 {
04283 if (logfile)
04284 {
04285
04286 FLOATWRITE(logfile, prob);
04287
04288 for (i=0; i<N; i++)
04289 FLOATWRITE(logfile, path_derivative_p(i,dim));
04290
04291 for (i=0; i<N; i++)
04292 FLOATWRITE(logfile, path_derivative_q(i,dim));
04293
04294 for (i=0; i<N; i++)
04295 for (j=0; j<N; j++)
04296 FLOATWRITE(logfile, path_derivative_a(i,j,dim));
04297
04298 for (i=0; i<N; i++)
04299 for (j=0; j<M; j++)
04300 FLOATWRITE(logfile, path_derivative_b(i,j,dim));
04301
04302 }
04303 }
04304 else
04305 {
04306 if (logfile)
04307 {
04308
04309 FLOATWRITE(logfile, prob);
04310
04311 for (i=0; model->get_learn_p(i)>=0; i++)
04312 FLOATWRITE(logfile, path_derivative_p(model->get_learn_p(i),dim));
04313
04314 for (i=0; model->get_learn_q(i)>=0; i++)
04315 FLOATWRITE(logfile, path_derivative_q(model->get_learn_q(i),dim));
04316
04317 for (q=0; model->get_learn_a(q,0)>=0; q++)
04318 {
04319 i=model->get_learn_a(q,0) ;
04320 j=model->get_learn_a(q,1) ;
04321 FLOATWRITE(logfile, path_derivative_a(i,j,dim));
04322 } ;
04323
04324 for (q=0; model->get_learn_b(q,0)>=0; q++)
04325 {
04326 i=model->get_learn_b(q,0) ;
04327 j=model->get_learn_b(q,1) ;
04328 FLOATWRITE(logfile, path_derivative_b(i,j,dim));
04329 } ;
04330 }
04331 } ;
04332 }
04333 save_model_bin(logfile) ;
04334
04335 result=true;
04336 SG_PRINT( "\n") ;
04337 return result;
04338 }
04339
04340 bool CHMM::save_model_derivatives_bin(FILE* file)
04341 {
04342 bool result=false;
04343 int32_t dim,i,j,q ;
04344 int32_t num_floats=0 ;
04345
04346 if (!model)
04347 SG_WARNING( "No definitions loaded -- writing derivatives of all weights\n") ;
04348 else
04349 SG_INFO( "writing derivatives of changed weights only\n") ;
04350
04351 #ifdef USE_HMMPARALLEL
04352 int32_t num_threads = parallel->get_num_threads();
04353 pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
04354 S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);
04355
04356 if (p_observations->get_num_vectors()<num_threads)
04357 num_threads=p_observations->get_num_vectors();
04358 #endif
04359
04360 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
04361 {
04362 if (dim%20==0)
04363 {
04364 SG_PRINT( ".") ;
04365
04366 } ;
04367
04368 #ifdef USE_HMMPARALLEL
04369 if (dim%num_threads==0)
04370 {
04371 for (i=0; i<num_threads; i++)
04372 {
04373 if (dim+i<p_observations->get_num_vectors())
04374 {
04375 params[i].hmm=this ;
04376 params[i].dim=dim+i ;
04377 pthread_create(&threads[i], NULL, bw_dim_prefetch, (void*)¶ms[i]) ;
04378 }
04379 }
04380
04381 for (i=0; i<num_threads; i++)
04382 {
04383 if (dim+i<p_observations->get_num_vectors())
04384 pthread_join(threads[i], NULL);
04385 }
04386 }
04387 #endif
04388
04389 float64_t prob=model_probability(dim) ;
04390 if (!model)
04391 {
04392 if (file)
04393 {
04394
04395 FLOATWRITE(file, prob);
04396
04397
04398 for (i=0; i<N; i++)
04399 FLOATWRITE(file, model_derivative_p(i,dim));
04400
04401 for (i=0; i<N; i++)
04402 FLOATWRITE(file, model_derivative_q(i,dim));
04403
04404
04405 for (i=0; i<N; i++)
04406 for (j=0; j<N; j++)
04407 FLOATWRITE(file, model_derivative_a(i,j,dim));
04408
04409 for (i=0; i<N; i++)
04410 for (j=0; j<M; j++)
04411 FLOATWRITE(file, model_derivative_b(i,j,dim));
04412
04413 if (dim==0)
04414 SG_INFO("Number of parameters (including posterior prob.): %i\n", num_floats) ;
04415 } ;
04416 }
04417 else
04418 {
04419 if (file)
04420 {
04421
04422 FLOATWRITE(file, prob);
04423
04424 for (i=0; model->get_learn_p(i)>=0; i++)
04425 FLOATWRITE(file, model_derivative_p(model->get_learn_p(i),dim));
04426
04427 for (i=0; model->get_learn_q(i)>=0; i++)
04428 FLOATWRITE(file, model_derivative_q(model->get_learn_q(i),dim));
04429
04430
04431 for (q=0; model->get_learn_a(q,1)>=0; q++)
04432 {
04433 i=model->get_learn_a(q,0) ;
04434 j=model->get_learn_a(q,1) ;
04435 FLOATWRITE(file, model_derivative_a(i,j,dim));
04436 } ;
04437
04438 for (q=0; model->get_learn_b(q,0)>=0; q++)
04439 {
04440 i=model->get_learn_b(q,0) ;
04441 j=model->get_learn_b(q,1) ;
04442 FLOATWRITE(file, model_derivative_b(i,j,dim));
04443 } ;
04444 if (dim==0)
04445 SG_INFO("Number of parameters (including posterior prob.): %i\n", num_floats) ;
04446 } ;
04447 } ;
04448 }
04449 save_model_bin(file) ;
04450
04451 #ifdef USE_HMMPARALLEL
04452 SG_FREE(threads);
04453 SG_FREE(params);
04454 #endif
04455
04456 result=true;
04457 SG_PRINT( "\n") ;
04458 return result;
04459 }
04460
04461 bool CHMM::save_model_derivatives(FILE* file)
04462 {
04463 bool result=false;
04464 int32_t dim,i,j;
04465
04466 if (file)
04467 {
04468
04469 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");
04470 fprintf(file,"%% dPr[...]=[ [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
04471 fprintf(file,"%% [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
04472 fprintf(file,"%% ............................. \n");
04473 fprintf(file,"%% [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_MM]\n");
04474 fprintf(file,"%% ];\n%%\n\nlog(dPr) = [\n");
04475
04476
04477 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
04478 {
04479 fprintf(file, "[ ");
04480
04481
04482 for (i=0; i<N; i++)
04483 fprintf(file,"%e, ", (double) model_derivative_p(i, dim) );
04484
04485 for (i=0; i<N; i++)
04486 fprintf(file,"%e, ", (double) model_derivative_q(i, dim) );
04487
04488
04489 for (i=0; i<N; i++)
04490 for (j=0; j<N; j++)
04491 fprintf(file, "%e,", (double) model_derivative_a(i,j,dim) );
04492
04493 for (i=0; i<N; i++)
04494 for (j=0; j<M; j++)
04495 fprintf(file, "%e,", (double) model_derivative_b(i,j,dim) );
04496
04497 fseek(file,ftell(file)-1,SEEK_SET);
04498 fprintf(file, " ];\n");
04499 }
04500
04501
04502 fprintf(file, "];");
04503
04504 result=true;
04505 }
04506 return result;
04507 }
04508
04509 bool CHMM::check_model_derivatives_combined()
04510 {
04511
04512 const float64_t delta=5e-4 ;
04513
04514 int32_t i ;
04515
04516
04517
04518
04519
04520
04521
04522
04523
04524
04525
04526
04527
04528
04529
04530
04531
04532
04533
04534
04535
04536
04537
04538
04539
04540
04541
04542
04543
04544
04545
04546 i=0;
04547 {
04548 for (int32_t j=0; j<M; j++)
04549 {
04550 float64_t old_b=get_b(i,j) ;
04551
04552 set_b(i,j, log(exp(old_b)-delta)) ;
04553 invalidate_model() ;
04554 float64_t prob_old=(model_probability(-1)*p_observations->get_num_vectors()) ;
04555
04556 set_b(i,j, log(exp(old_b)+delta)) ;
04557 invalidate_model() ;
04558 float64_t prob_new=(model_probability(-1)*p_observations->get_num_vectors());
04559
04560 float64_t deriv = (prob_new-prob_old)/(2*delta) ;
04561
04562 set_b(i,j, old_b) ;
04563 invalidate_model() ;
04564
04565 float64_t deriv_calc=0 ;
04566 for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
04567 {
04568 deriv_calc+=exp(model_derivative_b(i, j, dim)-model_probability(dim)) ;
04569 if (j==1)
04570 SG_INFO("deriv_calc[%i]=%e\n",dim,deriv_calc) ;
04571 } ;
04572
04573 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);
04574 } ;
04575 } ;
04576 return true ;
04577 }
04578
04579 bool CHMM::check_model_derivatives()
04580 {
04581 bool result=false;
04582 const float64_t delta=3e-4 ;
04583
04584 for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
04585 {
04586 int32_t i ;
04587
04588 for (i=0; i<N; i++)
04589 {
04590 for (int32_t j=0; j<N; j++)
04591 {
04592 float64_t old_a=get_a(i,j) ;
04593
04594 set_a(i,j, log(exp(old_a)-delta)) ;
04595 invalidate_model() ;
04596 float64_t prob_old=exp(model_probability(dim)) ;
04597
04598 set_a(i,j, log(exp(old_a)+delta)) ;
04599 invalidate_model() ;
04600 float64_t prob_new=exp(model_probability(dim));
04601
04602 float64_t deriv = (prob_new-prob_old)/(2*delta) ;
04603
04604 set_a(i,j, old_a) ;
04605 invalidate_model() ;
04606 float64_t deriv_calc=exp(model_derivative_a(i, j, dim)) ;
04607
04608 SG_DEBUG( "da(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc);
04609 invalidate_model() ;
04610 } ;
04611 } ;
04612 for (i=0; i<N; i++)
04613 {
04614 for (int32_t j=0; j<M; j++)
04615 {
04616 float64_t old_b=get_b(i,j) ;
04617
04618 set_b(i,j, log(exp(old_b)-delta)) ;
04619 invalidate_model() ;
04620 float64_t prob_old=exp(model_probability(dim)) ;
04621
04622 set_b(i,j, log(exp(old_b)+delta)) ;
04623 invalidate_model() ;
04624 float64_t prob_new=exp(model_probability(dim));
04625
04626 float64_t deriv = (prob_new-prob_old)/(2*delta) ;
04627
04628 set_b(i,j, old_b) ;
04629 invalidate_model() ;
04630 float64_t deriv_calc=exp(model_derivative_b(i, j, dim));
04631
04632 SG_DEBUG( "db(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/(deriv_calc));
04633 } ;
04634 } ;
04635
04636 #ifdef TEST
04637 for (i=0; i<N; i++)
04638 {
04639 float64_t old_p=get_p(i) ;
04640
04641 set_p(i, log(exp(old_p)-delta)) ;
04642 invalidate_model() ;
04643 float64_t prob_old=exp(model_probability(dim)) ;
04644
04645 set_p(i, log(exp(old_p)+delta)) ;
04646 invalidate_model() ;
04647 float64_t prob_new=exp(model_probability(dim));
04648 float64_t deriv = (prob_new-prob_old)/(2*delta) ;
04649
04650 set_p(i, old_p) ;
04651 invalidate_model() ;
04652 float64_t deriv_calc=exp(model_derivative_p(i, dim));
04653
04654
04655 SG_DEBUG( "dp(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc);
04656 } ;
04657 for (i=0; i<N; i++)
04658 {
04659 float64_t old_q=get_q(i) ;
04660
04661 set_q(i, log(exp(old_q)-delta)) ;
04662 invalidate_model() ;
04663 float64_t prob_old=exp(model_probability(dim)) ;
04664
04665 set_q(i, log(exp(old_q)+delta)) ;
04666 invalidate_model() ;
04667 float64_t prob_new=exp(model_probability(dim));
04668
04669 float64_t deriv = (prob_new-prob_old)/(2*delta) ;
04670
04671 set_q(i, old_q) ;
04672 invalidate_model() ;
04673 float64_t deriv_calc=exp(model_derivative_q(i, dim));
04674
04675
04676 SG_DEBUG( "dq(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc);
04677 } ;
04678 #endif
04679 }
04680 return result;
04681 }
04682
04683 #ifdef USE_HMMDEBUG
04684 bool CHMM::check_path_derivatives()
04685 {
04686 bool result=false;
04687 const float64_t delta=1e-4 ;
04688
04689 for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
04690 {
04691 int32_t i ;
04692
04693 for (i=0; i<N; i++)
04694 {
04695 for (int32_t j=0; j<N; j++)
04696 {
04697 float64_t old_a=get_a(i,j) ;
04698
04699 set_a(i,j, log(exp(old_a)-delta)) ;
04700 invalidate_model() ;
04701 float64_t prob_old=best_path(dim) ;
04702
04703 set_a(i,j, log(exp(old_a)+delta)) ;
04704 invalidate_model() ;
04705 float64_t prob_new=best_path(dim);
04706
04707 float64_t deriv = (prob_new-prob_old)/(2*delta) ;
04708
04709 set_a(i,j, old_a) ;
04710 invalidate_model() ;
04711 float64_t deriv_calc=path_derivative_a(i, j, dim) ;
04712
04713 SG_DEBUG( "da(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc);
04714 } ;
04715 } ;
04716 for (i=0; i<N; i++)
04717 {
04718 for (int32_t j=0; j<M; j++)
04719 {
04720 float64_t old_b=get_b(i,j) ;
04721
04722 set_b(i,j, log(exp(old_b)-delta)) ;
04723 invalidate_model() ;
04724 float64_t prob_old=best_path(dim) ;
04725
04726 set_b(i,j, log(exp(old_b)+delta)) ;
04727 invalidate_model() ;
04728 float64_t prob_new=best_path(dim);
04729
04730 float64_t deriv = (prob_new-prob_old)/(2*delta) ;
04731
04732 set_b(i,j, old_b) ;
04733 invalidate_model() ;
04734 float64_t deriv_calc=path_derivative_b(i, j, dim);
04735
04736 SG_DEBUG( "db(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/(deriv_calc));
04737 } ;
04738 } ;
04739
04740 for (i=0; i<N; i++)
04741 {
04742 float64_t old_p=get_p(i) ;
04743
04744 set_p(i, log(exp(old_p)-delta)) ;
04745 invalidate_model() ;
04746 float64_t prob_old=best_path(dim) ;
04747
04748 set_p(i, log(exp(old_p)+delta)) ;
04749 invalidate_model() ;
04750 float64_t prob_new=best_path(dim);
04751 float64_t deriv = (prob_new-prob_old)/(2*delta) ;
04752
04753 set_p(i, old_p) ;
04754 invalidate_model() ;
04755 float64_t deriv_calc=path_derivative_p(i, dim);
04756
04757
04758 SG_DEBUG( "dp(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc);
04759 } ;
04760 for (i=0; i<N; i++)
04761 {
04762 float64_t old_q=get_q(i) ;
04763
04764 set_q(i, log(exp(old_q)-delta)) ;
04765 invalidate_model() ;
04766 float64_t prob_old=best_path(dim) ;
04767
04768 set_q(i, log(exp(old_q)+delta)) ;
04769 invalidate_model() ;
04770 float64_t prob_new=best_path(dim);
04771
04772 float64_t deriv = (prob_new-prob_old)/(2*delta) ;
04773
04774 set_q(i, old_q) ;
04775 invalidate_model() ;
04776 float64_t deriv_calc=path_derivative_q(i, dim);
04777
04778
04779 SG_DEBUG( "dq(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc);
04780 } ;
04781 }
04782 return result;
04783 }
04784 #endif // USE_HMMDEBUG
04785
04786
04787 void CHMM::normalize(bool keep_dead_states)
04788 {
04789 int32_t i,j;
04790 const float64_t INF=-1e10;
04791 float64_t sum_p =INF;
04792
04793 for (i=0; i<N; i++)
04794 {
04795 sum_p=CMath::logarithmic_sum(sum_p, get_p(i));
04796
04797 float64_t sum_b =INF;
04798 float64_t sum_a =get_q(i);
04799
04800 for (j=0; j<N; j++)
04801 sum_a=CMath::logarithmic_sum(sum_a, get_a(i,j));
04802
04803 if (sum_a>CMath::ALMOST_NEG_INFTY/N || (!keep_dead_states) )
04804 {
04805 for (j=0; j<N; j++)
04806 set_a(i,j, get_a(i,j)-sum_a);
04807 set_q(i, get_q(i)-sum_a);
04808 }
04809
04810 for (j=0; j<M; j++)
04811 sum_b=CMath::logarithmic_sum(sum_b, get_b(i,j));
04812 for (j=0; j<M; j++)
04813 set_b(i,j, get_b(i,j)-sum_b);
04814 }
04815
04816 for (i=0; i<N; i++)
04817 set_p(i, get_p(i)-sum_p);
04818
04819 invalidate_model();
04820 }
04821
04822 bool CHMM::append_model(CHMM* app_model)
04823 {
04824 bool result=false;
04825 const int32_t num_states=app_model->get_N();
04826 int32_t i,j;
04827
04828 SG_DEBUG( "cur N:%d M:%d\n", N, M);
04829 SG_DEBUG( "old N:%d M:%d\n", app_model->get_N(), app_model->get_M());
04830 if (app_model->get_M() == get_M())
04831 {
04832 float64_t* n_p=SG_MALLOC(float64_t, N+num_states);
04833 float64_t* n_q=SG_MALLOC(float64_t, N+num_states);
04834 float64_t* n_a=SG_MALLOC(float64_t, (N+num_states)*(N+num_states));
04835
04836 float64_t* n_b=SG_MALLOC(float64_t, (N+num_states)*M);
04837
04838
04839 for (i=0; i<N+num_states; i++)
04840 {
04841 n_p[i]=-CMath::INFTY;
04842 n_q[i]=-CMath::INFTY;
04843
04844 for (j=0; j<N+num_states; j++)
04845 n_a[(N+num_states)*i+j]=-CMath::INFTY;
04846
04847 for (j=0; j<M; j++)
04848 n_b[M*i+j]=-CMath::INFTY;
04849 }
04850
04851
04852
04853
04854
04855
04856 for (i=0; i<N; i++)
04857 {
04858 n_p[i]=get_p(i);
04859
04860 for (j=0; j<N; j++)
04861 n_a[(N+num_states)*j+i]=get_a(i,j);
04862
04863 for (j=0; j<M; j++)
04864 {
04865 n_b[M*i+j]=get_b(i,j);
04866 }
04867 }
04868
04869
04870 for (i=0; i<app_model->get_N(); i++)
04871 {
04872 n_q[i+N]=app_model->get_q(i);
04873
04874 for (j=0; j<app_model->get_N(); j++)
04875 n_a[(N+num_states)*(j+N)+(i+N)]=app_model->get_a(i,j);
04876 for (j=0; j<app_model->get_M(); j++)
04877 n_b[M*(i+N)+j]=app_model->get_b(i,j);
04878 }
04879
04880
04881
04882 for (i=0; i<N; i++)
04883 {
04884 for (j=N; j<N+num_states; j++)
04885 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]);
04886 }
04887
04888 free_state_dependend_arrays();
04889 N+=num_states;
04890
04891 alloc_state_dependend_arrays();
04892
04893
04894 SG_FREE(initial_state_distribution_p);
04895 SG_FREE(end_state_distribution_q);
04896 SG_FREE(transition_matrix_a);
04897 SG_FREE(observation_matrix_b);
04898
04899 transition_matrix_a=n_a;
04900 observation_matrix_b=n_b;
04901 initial_state_distribution_p=n_p;
04902 end_state_distribution_q=n_q;
04903
04904 SG_WARNING( "not normalizing anymore, call normalize_hmm to make sure the hmm is valid!!\n");
04906 invalidate_model();
04907 }
04908 else
04909 SG_ERROR( "number of observations is different for append model, doing nothing!\n");
04910
04911 return result;
04912 }
04913
04914 bool CHMM::append_model(CHMM* app_model, float64_t* cur_out, float64_t* app_out)
04915 {
04916 bool result=false;
04917 const int32_t num_states=app_model->get_N()+2;
04918 int32_t i,j;
04919
04920 if (app_model->get_M() == get_M())
04921 {
04922 float64_t* n_p=SG_MALLOC(float64_t, N+num_states);
04923 float64_t* n_q=SG_MALLOC(float64_t, N+num_states);
04924 float64_t* n_a=SG_MALLOC(float64_t, (N+num_states)*(N+num_states));
04925
04926 float64_t* n_b=SG_MALLOC(float64_t, (N+num_states)*M);
04927
04928
04929 for (i=0; i<N+num_states; i++)
04930 {
04931 n_p[i]=-CMath::INFTY;
04932 n_q[i]=-CMath::INFTY;
04933
04934 for (j=0; j<N+num_states; j++)
04935 n_a[(N+num_states)*j+i]=-CMath::INFTY;
04936
04937 for (j=0; j<M; j++)
04938 n_b[M*i+j]=-CMath::INFTY;
04939 }
04940
04941
04942
04943
04944
04945
04946 for (i=0; i<N; i++)
04947 {
04948 n_p[i]=get_p(i);
04949
04950 for (j=0; j<N; j++)
04951 n_a[(N+num_states)*j+i]=get_a(i,j);
04952
04953 for (j=0; j<M; j++)
04954 {
04955 n_b[M*i+j]=get_b(i,j);
04956 }
04957 }
04958
04959
04960 for (i=0; i<app_model->get_N(); i++)
04961 {
04962 n_q[i+N+2]=app_model->get_q(i);
04963
04964 for (j=0; j<app_model->get_N(); j++)
04965 n_a[(N+num_states)*(j+N+2)+(i+N+2)]=app_model->get_a(i,j);
04966 for (j=0; j<app_model->get_M(); j++)
04967 n_b[M*(i+N+2)+j]=app_model->get_b(i,j);
04968 }
04969
04970
04971
04972
04973 for (i=0; i<M; i++)
04974 {
04975 n_b[M*N+i]=cur_out[i];
04976 n_b[M*(N+1)+i]=app_out[i];
04977 }
04978
04979
04980 for (i=0; i<N+num_states; i++)
04981 {
04982
04983 if (i==N+1)
04984 n_a[(N+num_states)*i + N]=0;
04985
04986
04987
04988 if (i<N)
04989 n_a[(N+num_states)*N+i]=get_q(i);
04990
04991
04992
04993 if (i>=N+2)
04994 n_a[(N+num_states)*i+(N+1)]=app_model->get_p(i-(N+2));
04995 }
04996
04997 free_state_dependend_arrays();
04998 N+=num_states;
04999
05000 alloc_state_dependend_arrays();
05001
05002
05003 SG_FREE(initial_state_distribution_p);
05004 SG_FREE(end_state_distribution_q);
05005 SG_FREE(transition_matrix_a);
05006 SG_FREE(observation_matrix_b);
05007
05008 transition_matrix_a=n_a;
05009 observation_matrix_b=n_b;
05010 initial_state_distribution_p=n_p;
05011 end_state_distribution_q=n_q;
05012
05013 SG_WARNING( "not normalizing anymore, call normalize_hmm to make sure the hmm is valid!!\n");
05015 invalidate_model();
05016 }
05017
05018 return result;
05019 }
05020
05021
05022 void CHMM::add_states(int32_t num_states, float64_t default_value)
05023 {
05024 int32_t i,j;
05025 const float64_t MIN_RAND=1e-2;
05026 const float64_t MAX_RAND=2e-1;
05027
05028 float64_t* n_p=SG_MALLOC(float64_t, N+num_states);
05029 float64_t* n_q=SG_MALLOC(float64_t, N+num_states);
05030 float64_t* n_a=SG_MALLOC(float64_t, (N+num_states)*(N+num_states));
05031
05032 float64_t* n_b=SG_MALLOC(float64_t, (N+num_states)*M);
05033
05034
05035
05036 for (i=0; i<N; i++)
05037 {
05038 n_p[i]=get_p(i);
05039 n_q[i]=get_q(i);
05040
05041 for (j=0; j<N; j++)
05042 n_a[(N+num_states)*j+i]=get_a(i,j);
05043
05044 for (j=0; j<M; j++)
05045 n_b[M*i+j]=get_b(i,j);
05046 }
05047
05048 for (i=N; i<N+num_states; i++)
05049 {
05050 n_p[i]=VAL_MACRO;
05051 n_q[i]=VAL_MACRO;
05052
05053 for (j=0; j<N; j++)
05054 n_a[(N+num_states)*i+j]=VAL_MACRO;
05055
05056 for (j=0; j<N+num_states; j++)
05057 n_a[(N+num_states)*j+i]=VAL_MACRO;
05058
05059 for (j=0; j<M; j++)
05060 n_b[M*i+j]=VAL_MACRO;
05061 }
05062 free_state_dependend_arrays();
05063 N+=num_states;
05064
05065 alloc_state_dependend_arrays();
05066
05067
05068 SG_FREE(initial_state_distribution_p);
05069 SG_FREE(end_state_distribution_q);
05070 SG_FREE(transition_matrix_a);
05071 SG_FREE(observation_matrix_b);
05072
05073 transition_matrix_a=n_a;
05074 observation_matrix_b=n_b;
05075 initial_state_distribution_p=n_p;
05076 end_state_distribution_q=n_q;
05077
05078 invalidate_model();
05079 normalize();
05080 }
05081
05082 void CHMM::chop(float64_t value)
05083 {
05084 for (int32_t i=0; i<N; i++)
05085 {
05086 int32_t j;
05087
05088 if (exp(get_p(i)) < value)
05089 set_p(i, CMath::ALMOST_NEG_INFTY);
05090
05091 if (exp(get_q(i)) < value)
05092 set_q(i, CMath::ALMOST_NEG_INFTY);
05093
05094 for (j=0; j<N; j++)
05095 {
05096 if (exp(get_a(i,j)) < value)
05097 set_a(i,j, CMath::ALMOST_NEG_INFTY);
05098 }
05099
05100 for (j=0; j<M; j++)
05101 {
05102 if (exp(get_b(i,j)) < value)
05103 set_b(i,j, CMath::ALMOST_NEG_INFTY);
05104 }
05105 }
05106 normalize();
05107 invalidate_model();
05108 }
05109
05110 bool CHMM::linear_train(bool right_align)
05111 {
05112 if (p_observations)
05113 {
05114 int32_t histsize=(get_M()*get_N());
05115 int32_t* hist=SG_MALLOC(int32_t, histsize);
05116 int32_t* startendhist=SG_MALLOC(int32_t, get_N());
05117 int32_t i,dim;
05118
05119 ASSERT(p_observations->get_max_vector_length()<=get_N());
05120
05121 for (i=0; i<histsize; i++)
05122 hist[i]=0;
05123
05124 for (i=0; i<get_N(); i++)
05125 startendhist[i]=0;
05126
05127 if (right_align)
05128 {
05129 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
05130 {
05131 int32_t len=0;
05132 bool free_vec;
05133 uint16_t* obs=p_observations->get_feature_vector(dim, len, free_vec);
05134
05135 ASSERT(len<=get_N());
05136 startendhist[(get_N()-len)]++;
05137
05138 for (i=0;i<len;i++)
05139 hist[(get_N()-len+i)*get_M() + *obs++]++;
05140
05141 p_observations->free_feature_vector(obs, dim, free_vec);
05142 }
05143
05144 set_q(get_N()-1, 1);
05145 for (i=0; i<get_N()-1; i++)
05146 set_q(i, 0);
05147
05148 for (i=0; i<get_N(); i++)
05149 set_p(i, startendhist[i]+PSEUDO);
05150
05151 for (i=0;i<get_N();i++)
05152 {
05153 for (int32_t j=0; j<get_N(); j++)
05154 {
05155 if (i==j-1)
05156 set_a(i,j, 1);
05157 else
05158 set_a(i,j, 0);
05159 }
05160 }
05161 }
05162 else
05163 {
05164 for (dim=0; dim<p_observations->get_num_vectors(); dim++)
05165 {
05166 int32_t len=0;
05167 bool free_vec;
05168 uint16_t* obs=p_observations->get_feature_vector(dim, len, free_vec);
05169
05170 ASSERT(len<=get_N());
05171 for (i=0;i<len;i++)
05172 hist[i*get_M() + *obs++]++;
05173
05174 startendhist[len-1]++;
05175
05176 p_observations->free_feature_vector(obs, dim, free_vec);
05177 }
05178
05179 set_p(0, 1);
05180 for (i=1; i<get_N(); i++)
05181 set_p(i, 0);
05182
05183 for (i=0; i<get_N(); i++)
05184 set_q(i, startendhist[i]+PSEUDO);
05185
05186 int32_t total=p_observations->get_num_vectors();
05187
05188 for (i=0;i<get_N();i++)
05189 {
05190 total-= startendhist[i] ;
05191
05192 for (int32_t j=0; j<get_N(); j++)
05193 {
05194 if (i==j-1)
05195 set_a(i,j, total+PSEUDO);
05196 else
05197 set_a(i,j, 0);
05198 }
05199 }
05200 ASSERT(total==0);
05201 }
05202
05203 for (i=0;i<get_N();i++)
05204 {
05205 for (int32_t j=0; j<get_M(); j++)
05206 {
05207 float64_t sum=0;
05208 int32_t offs=i*get_M()+ p_observations->get_masked_symbols((uint16_t) j, (uint8_t) 254);
05209
05210 for (int32_t k=0; k<p_observations->get_original_num_symbols(); k++)
05211 sum+=hist[offs+k];
05212
05213 set_b(i,j, (PSEUDO+hist[i*get_M()+j])/(sum+PSEUDO*p_observations->get_original_num_symbols()));
05214 }
05215 }
05216
05217 SG_FREE(hist);
05218 SG_FREE(startendhist);
05219 convert_to_log();
05220 invalidate_model();
05221 return true;
05222 }
05223 else
05224 return false;
05225 }
05226
05227 void CHMM::set_observation_nocache(CStringFeatures<uint16_t>* obs)
05228 {
05229 ASSERT(obs);
05230 p_observations=obs;
05231 SG_REF(obs);
05232
05233 if (obs)
05234 if (obs->get_num_symbols() > M)
05235 SG_ERROR( "number of symbols in observation (%ld) larger than M (%d)\n", (long) obs->get_num_symbols(), M);
05236
05237 if (!reused_caches)
05238 {
05239 #ifdef USE_HMMPARALLEL_STRUCTURES
05240 for (int32_t i=0; i<parallel->get_num_threads(); i++)
05241 {
05242 SG_FREE(alpha_cache[i].table);
05243 SG_FREE(beta_cache[i].table);
05244 SG_FREE(states_per_observation_psi[i]);
05245 SG_FREE(path[i]);
05246
05247 alpha_cache[i].table=NULL;
05248 beta_cache[i].table=NULL;
05249 states_per_observation_psi[i]=NULL;
05250 path[i]=NULL;
05251 } ;
05252 #else
05253 SG_FREE(alpha_cache.table);
05254 SG_FREE(beta_cache.table);
05255 SG_FREE(states_per_observation_psi);
05256 SG_FREE(path);
05257
05258 alpha_cache.table=NULL;
05259 beta_cache.table=NULL;
05260 states_per_observation_psi=NULL;
05261 path=NULL;
05262
05263 #endif //USE_HMMPARALLEL_STRUCTURES
05264 }
05265
05266 invalidate_model();
05267 }
05268
05269 void CHMM::set_observations(CStringFeatures<uint16_t>* obs, CHMM* lambda)
05270 {
05271 ASSERT(obs);
05272 SG_REF(obs);
05273 p_observations=obs;
05274
05275
05276
05277
05278 SG_REF(obs);
05279 features=obs;
05280
05281 SG_DEBUG("num symbols alphabet: %ld\n", obs->get_alphabet()->get_num_symbols());
05282 SG_DEBUG("num symbols: %ld\n", obs->get_num_symbols());
05283 SG_DEBUG("M: %d\n", M);
05284
05285 if (obs)
05286 {
05287 if (obs->get_num_symbols() > M)
05288 {
05289 SG_ERROR( "number of symbols in observation (%ld) larger than M (%d)\n", (long) obs->get_num_symbols(), M);
05290 }
05291 }
05292
05293 if (!reused_caches)
05294 {
05295 #ifdef USE_HMMPARALLEL_STRUCTURES
05296 for (int32_t i=0; i<parallel->get_num_threads(); i++)
05297 {
05298 SG_FREE(alpha_cache[i].table);
05299 SG_FREE(beta_cache[i].table);
05300 SG_FREE(states_per_observation_psi[i]);
05301 SG_FREE(path[i]);
05302
05303 alpha_cache[i].table=NULL;
05304 beta_cache[i].table=NULL;
05305 states_per_observation_psi[i]=NULL;
05306 path[i]=NULL;
05307 } ;
05308 #else
05309 SG_FREE(alpha_cache.table);
05310 SG_FREE(beta_cache.table);
05311 SG_FREE(states_per_observation_psi);
05312 SG_FREE(path);
05313
05314 alpha_cache.table=NULL;
05315 beta_cache.table=NULL;
05316 states_per_observation_psi=NULL;
05317 path=NULL;
05318
05319 #endif //USE_HMMPARALLEL_STRUCTURES
05320 }
05321
05322 if (obs!=NULL)
05323 {
05324 int32_t max_T=obs->get_max_vector_length();
05325
05326 if (lambda)
05327 {
05328 #ifdef USE_HMMPARALLEL_STRUCTURES
05329 for (int32_t i=0; i<parallel->get_num_threads(); i++)
05330 {
05331 this->alpha_cache[i].table= lambda->alpha_cache[i].table;
05332 this->beta_cache[i].table= lambda->beta_cache[i].table;
05333 this->states_per_observation_psi[i]=lambda->states_per_observation_psi[i] ;
05334 this->path[i]=lambda->path[i];
05335 } ;
05336 #else
05337 this->alpha_cache.table= lambda->alpha_cache.table;
05338 this->beta_cache.table= lambda->beta_cache.table;
05339 this->states_per_observation_psi= lambda->states_per_observation_psi;
05340 this->path=lambda->path;
05341 #endif //USE_HMMPARALLEL_STRUCTURES
05342
05343 this->reused_caches=true;
05344 }
05345 else
05346 {
05347 this->reused_caches=false;
05348 #ifdef USE_HMMPARALLEL_STRUCTURES
05349 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);
05350 for (int32_t i=0; i<parallel->get_num_threads(); i++)
05351 {
05352 if ((states_per_observation_psi[i]=SG_MALLOC(T_STATES,max_T*N))!=NULL)
05353 SG_DEBUG( "path_table[%i] successfully allocated\n",i) ;
05354 else
05355 SG_ERROR( "failed allocating memory for path_table[%i].\n",i) ;
05356 path[i]=SG_MALLOC(T_STATES, max_T);
05357 }
05358 #else // no USE_HMMPARALLEL_STRUCTURES
05359 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);
05360 if ((states_per_observation_psi=SG_MALLOC(T_STATES,max_T*N)) != NULL)
05361 SG_DONE();
05362 else
05363 SG_ERROR( "failed.\n") ;
05364
05365 path=SG_MALLOC(T_STATES, max_T);
05366 #endif // USE_HMMPARALLEL_STRUCTURES
05367 #ifdef USE_HMMCACHE
05368 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);
05369
05370 #ifdef USE_HMMPARALLEL_STRUCTURES
05371 for (int32_t i=0; i<parallel->get_num_threads(); i++)
05372 {
05373 if ((alpha_cache[i].table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N))!=NULL)
05374 SG_DEBUG( "alpha_cache[%i].table successfully allocated\n",i) ;
05375 else
05376 SG_ERROR("allocation of alpha_cache[%i].table failed\n",i) ;
05377
05378 if ((beta_cache[i].table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N)) != NULL)
05379 SG_DEBUG("beta_cache[%i].table successfully allocated\n",i) ;
05380 else
05381 SG_ERROR("allocation of beta_cache[%i].table failed\n",i) ;
05382 } ;
05383 #else // USE_HMMPARALLEL_STRUCTURES
05384 if ((alpha_cache.table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N)) != NULL)
05385 SG_DEBUG( "alpha_cache.table successfully allocated\n") ;
05386 else
05387 SG_ERROR( "allocation of alpha_cache.table failed\n") ;
05388
05389 if ((beta_cache.table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N)) != NULL)
05390 SG_DEBUG( "beta_cache.table successfully allocated\n") ;
05391 else
05392 SG_ERROR( "allocation of beta_cache.table failed\n") ;
05393
05394 #endif // USE_HMMPARALLEL_STRUCTURES
05395 #else // USE_HMMCACHE
05396 #ifdef USE_HMMPARALLEL_STRUCTURES
05397 for (int32_t i=0; i<parallel->get_num_threads(); i++)
05398 {
05399 alpha_cache[i].table=NULL ;
05400 beta_cache[i].table=NULL ;
05401 } ;
05402 #else //USE_HMMPARALLEL_STRUCTURES
05403 alpha_cache.table=NULL ;
05404 beta_cache.table=NULL ;
05405 #endif //USE_HMMPARALLEL_STRUCTURES
05406 #endif //USE_HMMCACHE
05407 }
05408 }
05409
05410
05411 invalidate_model();
05412 }
05413
05414 bool CHMM::permutation_entropy(int32_t window_width, int32_t sequence_number)
05415 {
05416 if (p_observations && window_width>0 &&
05417 ( sequence_number<0 || sequence_number < p_observations->get_num_vectors()))
05418 {
05419 int32_t min_sequence=sequence_number;
05420 int32_t max_sequence=sequence_number;
05421
05422 if (sequence_number<0)
05423 {
05424 min_sequence=0;
05425 max_sequence=p_observations->get_num_vectors();
05426 SG_INFO( "numseq: %d\n", max_sequence);
05427 }
05428
05429 SG_INFO( "min_sequence: %d max_sequence: %d\n", min_sequence, max_sequence);
05430 for (sequence_number=min_sequence; sequence_number<max_sequence; sequence_number++)
05431 {
05432 int32_t sequence_length=0;
05433 bool free_vec;
05434 uint16_t* obs=p_observations->get_feature_vector(sequence_number, sequence_length, free_vec);
05435
05436 int32_t histsize=get_M();
05437 int64_t* hist=SG_MALLOC(int64_t, histsize);
05438 int32_t i,j;
05439
05440 for (i=0; i<sequence_length-window_width; i++)
05441 {
05442 for (j=0; j<histsize; j++)
05443 hist[j]=0;
05444
05445 uint16_t* ptr=&obs[i];
05446 for (j=0; j<window_width; j++)
05447 {
05448 hist[*ptr++]++;
05449 }
05450
05451 float64_t perm_entropy=0;
05452 for (j=0; j<get_M(); j++)
05453 {
05454 float64_t p=
05455 (((float64_t) hist[j])+PSEUDO)/
05456 (window_width+get_M()*PSEUDO);
05457 perm_entropy+=p*log(p);
05458 }
05459
05460 SG_PRINT( "%f\n", perm_entropy);
05461 }
05462 p_observations->free_feature_vector(obs, sequence_number, free_vec);
05463
05464 SG_FREE(hist);
05465 }
05466 return true;
05467 }
05468 else
05469 return false;
05470 }
05471
05472 float64_t CHMM::get_log_derivative(int32_t num_param, int32_t num_example)
05473 {
05474 if (num_param<N)
05475 return model_derivative_p(num_param, num_example);
05476 else if (num_param<2*N)
05477 return model_derivative_q(num_param-N, num_example);
05478 else if (num_param<N*(N+2))
05479 {
05480 int32_t k=num_param-2*N;
05481 int32_t i=k/N;
05482 int32_t j=k%N;
05483 return model_derivative_a(i,j, num_example);
05484 }
05485 else if (num_param<N*(N+2+M))
05486 {
05487 int32_t k=num_param-N*(N+2);
05488 int32_t i=k/M;
05489 int32_t j=k%M;
05490 return model_derivative_b(i,j, num_example);
05491 }
05492
05493 ASSERT(false);
05494 return -1;
05495 }
05496
05497 float64_t CHMM::get_log_model_parameter(int32_t num_param)
05498 {
05499 if (num_param<N)
05500 return get_p(num_param);
05501 else if (num_param<2*N)
05502 return get_q(num_param-N);
05503 else if (num_param<N*(N+2))
05504 return transition_matrix_a[num_param-2*N];
05505 else if (num_param<N*(N+2+M))
05506 return observation_matrix_b[num_param-N*(N+2)];
05507
05508 ASSERT(false);
05509 return -1;
05510 }
05511
05512
05513
05514 bool CHMM::converged(float64_t x, float64_t y)
05515 {
05516 float64_t diff=y-x;
05517 float64_t absdiff=fabs(diff);
05518
05519 SG_INFO( "\n #%03d\tbest result so far: %G (eps: %f)", iteration_count, y, diff);
05520
05521 if (iteration_count--==0 || (absdiff<epsilon && conv_it<=0))
05522 {
05523 iteration_count=iterations;
05524 SG_INFO( "...finished\n");
05525 conv_it=5;
05526 return true;
05527 }
05528 else
05529 {
05530 if (absdiff<epsilon)
05531 conv_it--;
05532 else
05533 conv_it=5;
05534
05535 return false;
05536 }
05537 }
05538
05539 bool CHMM::baum_welch_viterbi_train(BaumWelchViterbiType type)
05540 {
05541 CHMM* estimate=new CHMM(this);
05542 CHMM* working=this;
05543 float64_t prob_max=-CMath::INFTY;
05544 float64_t prob=-CMath::INFTY;
05545 float64_t prob_train=CMath::ALMOST_NEG_INFTY;
05546 iteration_count=iterations;
05547
05548 while (!converged(prob, prob_train) && (!CSignal::cancel_computations()))
05549 {
05550 CMath::swap(working, estimate);
05551 prob=prob_train;
05552
05553 switch (type) {
05554 case BW_NORMAL:
05555 working->estimate_model_baum_welch(estimate); break;
05556 case BW_TRANS:
05557 working->estimate_model_baum_welch_trans(estimate); break;
05558 case BW_DEFINED:
05559 working->estimate_model_baum_welch_defined(estimate); break;
05560 case VIT_NORMAL:
05561 working->estimate_model_viterbi(estimate); break;
05562 case VIT_DEFINED:
05563 working->estimate_model_viterbi_defined(estimate); break;
05564 }
05565 prob_train=estimate->model_probability();
05566
05567 if (prob_max<prob_train)
05568 prob_max=prob_train;
05569 }
05570
05571 if (estimate == this)
05572 {
05573 estimate->copy_model(working);
05574 delete working;
05575 }
05576 else
05577 delete estimate;
05578
05579 return true;
05580 }