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