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

SHOGUN Machine Learning Toolbox - Documentation