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