HMM.h

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 
00012 #ifndef __CHMM_H__
00013 #define __CHMM_H__
00014 
00015 #include "lib/Mathematics.h"
00016 #include "lib/common.h"
00017 #include "lib/io.h"
00018 #include "lib/config.h"
00019 #include "features/Features.h"
00020 #include "features/StringFeatures.h"
00021 #include "distributions/Distribution.h"
00022 
00023 #include <stdio.h>
00024 
00025 #ifdef USE_HMMPARALLEL
00026 #define USE_HMMPARALLEL_STRUCTURES 1
00027 #endif
00028 
00029 namespace shogun
00030 {
00031     class CFeatures;
00032     template <class ST> class CStringFeatures;
00035 
00037 typedef float64_t T_ALPHA_BETA_TABLE;
00038 
00039 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00040 
00041 struct T_ALPHA_BETA
00042 {
00044     int32_t dimension;
00045 
00047     T_ALPHA_BETA_TABLE* table;
00048 
00050     bool updated;
00051 
00053     float64_t sum;
00054 };
00055 #endif // DOXYGEN_SHOULD_SKIP_THIS
00056 
00061 #ifdef USE_BIGSTATES
00062 typedef uint16_t T_STATES ;
00063 #else
00064 typedef uint8_t T_STATES ;
00065 #endif
00066 typedef T_STATES* P_STATES ;
00067 
00069 
00070 enum BaumWelchViterbiType
00071 {
00072     BW_NORMAL,
00073     BW_TRANS,
00074     BW_DEFINED,
00075     VIT_NORMAL,
00076     VIT_DEFINED
00077 };
00078 
00079 
00081 class CModel
00082 {
00083     public:
00085         CModel();
00086 
00088         virtual ~CModel();
00089 
00091         inline void sort_learn_a()
00092         {
00093             CMath::sort(learn_a,2) ;
00094         }
00095 
00097         inline void sort_learn_b()
00098         {
00099             CMath::sort(learn_b,2) ;
00100         }
00101 
00106 
00107         inline int32_t get_learn_a(int32_t line, int32_t column) const
00108         {
00109             return learn_a[line*2 + column];
00110         }
00111 
00113         inline int32_t get_learn_b(int32_t line, int32_t column) const 
00114         {
00115             return learn_b[line*2 + column];
00116         }
00117 
00119         inline int32_t get_learn_p(int32_t offset) const 
00120         {
00121             return learn_p[offset];
00122         }
00123 
00125         inline int32_t get_learn_q(int32_t offset) const 
00126         {
00127             return learn_q[offset];
00128         }
00129 
00131         inline int32_t get_const_a(int32_t line, int32_t column) const
00132         {
00133             return const_a[line*2 + column];
00134         }
00135 
00137         inline int32_t get_const_b(int32_t line, int32_t column) const 
00138         {
00139             return const_b[line*2 + column];
00140         }
00141 
00143         inline int32_t get_const_p(int32_t offset) const 
00144         {
00145             return const_p[offset];
00146         }
00147 
00149         inline int32_t get_const_q(int32_t offset) const
00150         {
00151             return const_q[offset];
00152         }
00153 
00155         inline float64_t get_const_a_val(int32_t line) const
00156         {
00157             return const_a_val[line];
00158         }
00159 
00161         inline float64_t get_const_b_val(int32_t line) const 
00162         {
00163             return const_b_val[line];
00164         }
00165 
00167         inline float64_t get_const_p_val(int32_t offset) const 
00168         {
00169             return const_p_val[offset];
00170         }
00171 
00173         inline float64_t get_const_q_val(int32_t offset) const
00174         {
00175             return const_q_val[offset];
00176         }
00177 #ifdef FIX_POS
00178 
00179         inline char get_fix_pos_state(int32_t pos, T_STATES state, T_STATES num_states)
00180         {
00181 #ifdef HMM_DEBUG
00182             if ((pos<0)||(pos*num_states+state>65336))
00183                 SG_DEBUG("index out of range in get_fix_pos_state(%i,%i,%i) \n", pos,state,num_states) ;
00184 #endif
00185             return fix_pos_state[pos*num_states+state] ;
00186         }
00187 #endif
00188 
00189 
00194 
00195         inline void set_learn_a(int32_t offset, int32_t value)
00196         {
00197             learn_a[offset]=value;
00198         }
00199 
00201         inline void set_learn_b(int32_t offset, int32_t value)
00202         {
00203             learn_b[offset]=value;
00204         }
00205 
00207         inline void set_learn_p(int32_t offset, int32_t value)
00208         {
00209             learn_p[offset]=value;
00210         }
00211 
00213         inline void set_learn_q(int32_t offset, int32_t value)
00214         {
00215             learn_q[offset]=value;
00216         }
00217 
00219         inline void set_const_a(int32_t offset, int32_t value)
00220         {
00221             const_a[offset]=value;
00222         }
00223 
00225         inline void set_const_b(int32_t offset, int32_t value)
00226         {
00227             const_b[offset]=value;
00228         }
00229 
00231         inline void set_const_p(int32_t offset, int32_t value)
00232         {
00233             const_p[offset]=value;
00234         }
00235 
00237         inline void set_const_q(int32_t offset, int32_t value)
00238         {
00239             const_q[offset]=value;
00240         }
00241 
00243         inline void set_const_a_val(int32_t offset, float64_t value)
00244         {
00245             const_a_val[offset]=value;
00246         }
00247 
00249         inline void set_const_b_val(int32_t offset, float64_t value)
00250         {
00251             const_b_val[offset]=value;
00252         }
00253 
00255         inline void set_const_p_val(int32_t offset, float64_t value)
00256         {
00257             const_p_val[offset]=value;
00258         }
00259 
00261         inline void set_const_q_val(int32_t offset, float64_t value)
00262         {
00263             const_q_val[offset]=value;
00264         }
00265 #ifdef FIX_POS
00266 
00267         inline void set_fix_pos_state(
00268             int32_t pos, T_STATES state, T_STATES num_states, char value)
00269         {
00270 #ifdef HMM_DEBUG
00271             if ((pos<0)||(pos*num_states+state>65336))
00272                 SG_DEBUG("index out of range in set_fix_pos_state(%i,%i,%i,%i) [%i]\n", pos,state,num_states,(int)value, pos*num_states+state) ;
00273 #endif
00274             fix_pos_state[pos*num_states+state]=value;
00275             if (value==FIX_ALLOWED)
00276                 for (int32_t i=0; i<num_states; i++)
00277                     if (get_fix_pos_state(pos,i,num_states)==FIX_DEFAULT)
00278                         set_fix_pos_state(pos,i,num_states,FIX_DISALLOWED) ;
00279         }
00281 
00283         const static char FIX_DISALLOWED ;
00284 
00286         const static char FIX_ALLOWED ;
00287 
00289         const static char FIX_DEFAULT ;
00290 
00292         const static float64_t DISALLOWED_PENALTY ;
00293 #endif
00294     protected:
00301 
00302         int32_t* learn_a;
00303 
00305         int32_t* learn_b;
00306 
00308         int32_t* learn_p;
00309 
00311         int32_t* learn_q;
00313 
00320 
00321         int32_t* const_a;
00322 
00324         int32_t* const_b;
00325 
00327         int32_t* const_p;
00328 
00330         int32_t* const_q;       
00331 
00332 
00334         float64_t* const_a_val;
00335 
00337         float64_t* const_b_val;
00338 
00340         float64_t* const_p_val;
00341 
00343         float64_t* const_q_val;     
00344 
00345 #ifdef FIX_POS
00346 
00349         char* fix_pos_state;
00350 #endif
00351 
00352 };
00353 
00354 
00365 class CHMM : public CDistribution
00366 {
00367     private:
00368 
00369         T_STATES trans_list_len ;
00370         T_STATES **trans_list_forward  ;
00371         T_STATES *trans_list_forward_cnt  ;
00372         float64_t **trans_list_forward_val ;
00373         T_STATES **trans_list_backward  ;
00374         T_STATES *trans_list_backward_cnt  ;
00375         bool mem_initialized ;
00376 
00377 #ifdef USE_HMMPARALLEL_STRUCTURES
00378 
00380         struct S_DIM_THREAD_PARAM
00381         {
00382             CHMM* hmm;
00383             int32_t dim;
00384             float64_t prob_sum;
00385         };
00386 
00388         struct S_BW_THREAD_PARAM
00389         {
00390             CHMM* hmm;
00391             int32_t dim_start;
00392             int32_t dim_stop;
00393 
00394             float64_t ret;
00395 
00396             float64_t* p_buf;
00397             float64_t* q_buf;
00398             float64_t* a_buf;
00399             float64_t* b_buf;
00400         };
00401 
00402         inline T_ALPHA_BETA & ALPHA_CACHE(int32_t dim) {
00403             return alpha_cache[dim%parallel->get_num_threads()] ; } ;
00404         inline T_ALPHA_BETA & BETA_CACHE(int32_t dim) {
00405             return beta_cache[dim%parallel->get_num_threads()] ; } ;
00406 #ifdef USE_LOGSUMARRAY 
00407         inline float64_t* ARRAYS(int32_t dim) {
00408             return arrayS[dim%parallel->get_num_threads()] ; } ;
00409 #endif
00410         inline float64_t* ARRAYN1(int32_t dim) {
00411             return arrayN1[dim%parallel->get_num_threads()] ; } ;
00412         inline float64_t* ARRAYN2(int32_t dim) {
00413             return arrayN2[dim%parallel->get_num_threads()] ; } ;
00414         inline T_STATES* STATES_PER_OBSERVATION_PSI(int32_t dim) {
00415             return states_per_observation_psi[dim%parallel->get_num_threads()] ; } ;
00416         inline const T_STATES* STATES_PER_OBSERVATION_PSI(int32_t dim) const {
00417             return states_per_observation_psi[dim%parallel->get_num_threads()] ; } ;
00418         inline T_STATES* PATH(int32_t dim) {
00419             return path[dim%parallel->get_num_threads()] ; } ;
00420         inline bool & PATH_PROB_UPDATED(int32_t dim) {
00421             return path_prob_updated[dim%parallel->get_num_threads()] ; } ;
00422         inline int32_t & PATH_PROB_DIMENSION(int32_t dim) {
00423             return path_prob_dimension[dim%parallel->get_num_threads()] ; } ;
00424 #else
00425         inline T_ALPHA_BETA & ALPHA_CACHE(int32_t /*dim*/) {
00426             return alpha_cache ; } ;
00427         inline T_ALPHA_BETA & BETA_CACHE(int32_t /*dim*/) {
00428             return beta_cache ; } ;
00429 #ifdef USE_LOGSUMARRAY
00430         inline float64_t* ARRAYS(int32_t dim) {
00431             return arrayS ; } ;
00432 #endif
00433         inline float64_t* ARRAYN1(int32_t /*dim*/) {
00434             return arrayN1 ; } ;
00435         inline float64_t* ARRAYN2(int32_t /*dim*/) {
00436             return arrayN2 ; } ;
00437         inline T_STATES* STATES_PER_OBSERVATION_PSI(int32_t /*dim*/) {
00438             return states_per_observation_psi ; } ;
00439         inline const T_STATES* STATES_PER_OBSERVATION_PSI(int32_t /*dim*/) const {
00440             return states_per_observation_psi ; } ;
00441         inline T_STATES* PATH(int32_t /*dim*/) {
00442             return path ; } ;
00443         inline bool & PATH_PROB_UPDATED(int32_t /*dim*/) {
00444             return path_prob_updated ; } ;
00445         inline int32_t & PATH_PROB_DIMENSION(int32_t /*dim*/) {
00446             return path_prob_dimension ; } ;
00447 #endif
00448 
00453         bool converged(float64_t x, float64_t y);
00454 
00460     public:
00471         CHMM(
00472             int32_t N, int32_t M, CModel* model, float64_t PSEUDO);
00473         CHMM(
00474             CStringFeatures<uint16_t>* obs, int32_t N, int32_t M,
00475             float64_t PSEUDO);
00476         CHMM(
00477             int32_t N, float64_t* p, float64_t* q, float64_t* a);
00478         CHMM(
00479             int32_t N, float64_t* p, float64_t* q, int32_t num_trans,
00480             float64_t* a_trans);
00481 
00486         CHMM(FILE* model_file, float64_t PSEUDO);
00487 
00489         CHMM(CHMM* h);
00490 
00492         virtual ~CHMM();
00493 
00502         virtual bool train(CFeatures* data=NULL);
00503         virtual inline int32_t get_num_model_parameters() { return N*(N+M+2); }
00504         virtual float64_t get_log_model_parameter(int32_t num_param);
00505         virtual float64_t get_log_derivative(int32_t num_param, int32_t num_example);
00506         virtual float64_t get_log_likelihood_example(int32_t num_example)
00507         {
00508             return model_probability(num_example);
00509         }
00510 
00516         bool initialize(CModel* model, float64_t PSEUDO, FILE* model_file=NULL);
00518 
00520         bool alloc_state_dependend_arrays();
00521 
00523         void free_state_dependend_arrays();
00524 
00536         float64_t forward_comp(int32_t time, int32_t state, int32_t dimension);
00537         float64_t forward_comp_old(
00538             int32_t time, int32_t state, int32_t dimension);
00539 
00547         float64_t backward_comp(int32_t time, int32_t state, int32_t dimension);
00548         float64_t backward_comp_old(
00549             int32_t time, int32_t state, int32_t dimension);
00550 
00555         float64_t best_path(int32_t dimension);
00556         inline uint16_t get_best_path_state(int32_t dim, int32_t t)
00557         {
00558             ASSERT(PATH(dim));
00559             return PATH(dim)[t];
00560         }
00561 
00564         float64_t model_probability_comp() ;
00565 
00567         inline float64_t model_probability(int32_t dimension=-1)
00568         {
00569             //for faster calculation cache model probability
00570             if (dimension==-1)
00571             {
00572                 if (mod_prob_updated)
00573                     return mod_prob/p_observations->get_num_vectors();
00574                 else
00575                     return model_probability_comp()/p_observations->get_num_vectors();
00576             }
00577             else
00578                 return forward(p_observations->get_vector_length(dimension), 0, dimension);
00579         }
00580 
00586         inline float64_t linear_model_probability(int32_t dimension)
00587         {
00588             float64_t lik=0;
00589             int32_t len=0;
00590             bool free_vec;
00591             uint16_t* o=p_observations->get_feature_vector(dimension, len, free_vec);
00592             float64_t* obs_b=observation_matrix_b;
00593 
00594             ASSERT(N==len);
00595 
00596             for (int32_t i=0; i<N; i++)
00597             {
00598                 lik+=obs_b[*o++];
00599                 obs_b+=M;
00600             }
00601             p_observations->free_feature_vector(o, dimension, free_vec);
00602             return lik;
00603 
00604             // sorry, the above code is the speed optimized version of :
00605             /*  float64_t lik=0;
00606 
00607                 for (int32_t i=0; i<N; i++)
00608                 lik+=get_b(i, p_observations->get_feature(dimension, i));
00609                 return lik;
00610                 */
00611             // : that
00612         }
00613 
00615 
00618         inline bool set_iterations(int32_t num) { iterations=num; return true; }
00619         inline int32_t get_iterations() { return iterations; }
00620         inline bool set_epsilon (float64_t eps) { epsilon=eps; return true; }
00621         inline float64_t get_epsilon() { return epsilon; }
00622 
00626         bool baum_welch_viterbi_train(BaumWelchViterbiType type);
00627 
00634         void estimate_model_baum_welch(CHMM* train);
00635         void estimate_model_baum_welch_trans(CHMM* train);
00636 
00637 #ifdef USE_HMMPARALLEL_STRUCTURES
00638         void ab_buf_comp(
00639             float64_t* p_buf, float64_t* q_buf, float64_t* a_buf,
00640             float64_t* b_buf, int32_t dim) ;
00641 #else
00642         void estimate_model_baum_welch_old(CHMM* train);
00643 #endif
00644 
00648         void estimate_model_baum_welch_defined(CHMM* train);
00649 
00653         void estimate_model_viterbi(CHMM* train);
00654 
00658         void estimate_model_viterbi_defined(CHMM* train);
00659 
00661 
00663         bool linear_train(bool right_align=false);
00664 
00666         bool permutation_entropy(int32_t window_width, int32_t sequence_number);
00667 
00674         void output_model(bool verbose=false);
00675 
00677         void output_model_defined(bool verbose=false);
00679 
00680 
00683 
00685         void normalize(bool keep_dead_states=false);
00686 
00690         void add_states(int32_t num_states, float64_t default_val=0);
00691 
00697         bool append_model(
00698             CHMM* append_model, float64_t* cur_out, float64_t* app_out);
00699 
00703         bool append_model(CHMM* append_model);
00704 
00706         void chop(float64_t value);
00707 
00709         void convert_to_log();
00710 
00712         void init_model_random();
00713 
00719         void init_model_defined();
00720 
00722         void clear_model();
00723 
00725         void clear_model_defined();
00726 
00728         void copy_model(CHMM* l);
00729 
00734         void invalidate_model();
00735 
00739         inline bool get_status() const 
00740         {   
00741             return status; 
00742         } 
00743 
00745         inline float64_t get_pseudo() const
00746         {
00747             return PSEUDO ;
00748         }
00749 
00751         inline void set_pseudo(float64_t pseudo) 
00752         {
00753             PSEUDO=pseudo ;
00754         }
00755 
00756 #ifdef USE_HMMPARALLEL_STRUCTURES
00757         static void* bw_dim_prefetch(void * params);
00758         static void* bw_single_dim_prefetch(void * params);
00759         static void* vit_dim_prefetch(void * params);
00760 #endif
00761 
00762 #ifdef FIX_POS
00763 
00766         inline bool set_fix_pos_state(int32_t pos, T_STATES state, char value)
00767         {
00768             if (!model)
00769                 return false ;
00770             model->set_fix_pos_state(pos, state, N, value) ;
00771             return true ;
00772         } ;
00773 #endif  
00774 
00775 
00784         void set_observations(CStringFeatures<uint16_t>* obs, CHMM* hmm=NULL);
00785 
00789         void set_observation_nocache(CStringFeatures<uint16_t>* obs);
00790 
00792         inline CStringFeatures<uint16_t>* get_observations()
00793         {
00794             SG_REF(p_observations);
00795             return p_observations;
00796         }
00798 
00866         bool load_definitions(FILE* file, bool verbose, bool initialize=true);
00867 
00903         bool load_model(FILE* file);
00904 
00908         bool save_model(FILE* file);
00909 
00913         bool save_model_derivatives(FILE* file);
00914 
00918         bool save_model_derivatives_bin(FILE* file);
00919 
00923         bool save_model_bin(FILE* file);
00924 
00926         bool check_model_derivatives() ;
00927         bool check_model_derivatives_combined() ;
00928 
00934         T_STATES* get_path(int32_t dim, float64_t& prob);
00935 
00939         bool save_path(FILE* file);
00940 
00944         bool save_path_derivatives(FILE* file);
00945 
00949         bool save_path_derivatives_bin(FILE* file);
00950 
00951 #ifdef USE_HMMDEBUG
00952 
00953         bool check_path_derivatives() ;
00954 #endif //USE_HMMDEBUG
00955 
00959         bool save_likelihood_bin(FILE* file);
00960 
00964         bool save_likelihood(FILE* file);
00966 
00972 
00974         inline T_STATES get_N() const { return N ; }
00975 
00977         inline int32_t get_M() const { return M ; }
00978 
00983         inline void set_q(T_STATES offset, float64_t value)
00984         {
00985 #ifdef HMM_DEBUG
00986             if (offset>=N)
00987                 SG_DEBUG("index out of range in set_q(%i,%e) [%i]\n", offset,value,N) ;
00988 #endif
00989             end_state_distribution_q[offset]=value;
00990         }
00991 
00996         inline void set_p(T_STATES offset, float64_t value)
00997         {
00998 #ifdef HMM_DEBUG
00999             if (offset>=N)
01000                 SG_DEBUG("index out of range in set_p(%i,.) [%i]\n", offset,N) ;
01001 #endif
01002             initial_state_distribution_p[offset]=value;
01003         }
01004 
01010         inline void set_A(T_STATES line_, T_STATES column, float64_t value)
01011         {
01012 #ifdef HMM_DEBUG
01013             if ((line_>N)||(column>N))
01014                 SG_DEBUG("index out of range in set_A(%i,%i,.) [%i,%i]\n",line_,column,N,N) ;
01015 #endif
01016             transition_matrix_A[line_+column*N]=value;
01017         }
01018 
01024         inline void set_a(T_STATES line_, T_STATES column, float64_t value)
01025         {
01026 #ifdef HMM_DEBUG
01027             if ((line_>N)||(column>N))
01028                 SG_DEBUG("index out of range in set_a(%i,%i,.) [%i,%i]\n",line_,column,N,N) ;
01029 #endif
01030             transition_matrix_a[line_+column*N]=value; // look also best_path!
01031         }
01032 
01038         inline void set_B(T_STATES line_, uint16_t column, float64_t value)
01039         {
01040 #ifdef HMM_DEBUG
01041             if ((line_>=N)||(column>=M))
01042                 SG_DEBUG("index out of range in set_B(%i,%i) [%i,%i]\n", line_, column,N,M) ;
01043 #endif
01044             observation_matrix_B[line_*M+column]=value;
01045         }
01046 
01052         inline void set_b(T_STATES line_, uint16_t column, float64_t value)
01053         {
01054 #ifdef HMM_DEBUG
01055             if ((line_>=N)||(column>=M))
01056                 SG_DEBUG("index out of range in set_b(%i,%i) [%i,%i]\n", line_, column,N,M) ;
01057 #endif
01058             observation_matrix_b[line_*M+column]=value;
01059         }
01060 
01067         inline void set_psi(
01068             int32_t time, T_STATES state, T_STATES value, int32_t dimension)
01069         {
01070 #ifdef HMM_DEBUG
01071             if ((time>=p_observations->get_max_vector_length())||(state>N))
01072                 SG_DEBUG("index out of range in set_psi(%i,%i,.) [%i,%i]\n",time,state,p_observations->get_max_vector_length(),N) ;
01073 #endif
01074             STATES_PER_OBSERVATION_PSI(dimension)[time*N+state]=value;
01075         }
01076 
01081         inline float64_t get_q(T_STATES offset) const 
01082         {
01083 #ifdef HMM_DEBUG
01084             if (offset>=N)
01085                 SG_DEBUG("index out of range in %e=get_q(%i) [%i]\n", end_state_distribution_q[offset],offset,N) ;
01086 #endif
01087             return end_state_distribution_q[offset];
01088         }
01089 
01094         inline float64_t get_p(T_STATES offset) const 
01095         {
01096 #ifdef HMM_DEBUG
01097             if (offset>=N)
01098                 SG_DEBUG("index out of range in get_p(%i,.) [%i]\n", offset,N) ;
01099 #endif
01100             return initial_state_distribution_p[offset];
01101         }
01102 
01108         inline float64_t get_A(T_STATES line_, T_STATES column) const
01109         {
01110 #ifdef HMM_DEBUG
01111             if ((line_>N)||(column>N))
01112                 SG_DEBUG("index out of range in get_A(%i,%i) [%i,%i]\n",line_,column,N,N) ;
01113 #endif
01114             return transition_matrix_A[line_+column*N];
01115         }
01116 
01122         inline float64_t get_a(T_STATES line_, T_STATES column) const
01123         {
01124 #ifdef HMM_DEBUG
01125             if ((line_>N)||(column>N))
01126                 SG_DEBUG("index out of range in get_a(%i,%i) [%i,%i]\n",line_,column,N,N) ;
01127 #endif
01128             return transition_matrix_a[line_+column*N]; // look also best_path()!
01129         }
01130 
01136         inline float64_t get_B(T_STATES line_, uint16_t column) const
01137         {
01138 #ifdef HMM_DEBUG
01139             if ((line_>=N)||(column>=M))
01140                 SG_DEBUG("index out of range in get_B(%i,%i) [%i,%i]\n", line_, column,N,M) ;
01141 #endif
01142             return observation_matrix_B[line_*M+column];
01143         }
01144 
01150         inline float64_t get_b(T_STATES line_, uint16_t column) const 
01151         {
01152 #ifdef HMM_DEBUG
01153             if ((line_>=N)||(column>=M))
01154                 SG_DEBUG("index out of range in get_b(%i,%i) [%i,%i]\n", line_, column,N,M) ;
01155 #endif
01156             //SG_PRINT("idx %d\n", line_*M+column);
01157             return observation_matrix_b[line_*M+column];
01158         }
01159 
01166         inline T_STATES get_psi(
01167             int32_t time, T_STATES state, int32_t dimension) const
01168         {
01169 #ifdef HMM_DEBUG
01170             if ((time>=p_observations->get_max_vector_length())||(state>N))
01171                 SG_DEBUG("index out of range in get_psi(%i,%i) [%i,%i]\n",time,state,p_observations->get_max_vector_length(),N) ;
01172 #endif
01173             return STATES_PER_OBSERVATION_PSI(dimension)[time*N+state];
01174         }
01175 
01177 
01179         inline virtual const char* get_name() const { return "HMM"; }
01180     protected:
01185 
01186         int32_t M;
01187 
01189         int32_t N;
01190 
01192         float64_t PSEUDO;
01193 
01194         // line number during processing input files
01195         int32_t line;
01196 
01198         CStringFeatures<uint16_t>* p_observations;
01199 
01200         //train definition for HMM
01201         CModel* model;
01202 
01204         float64_t* transition_matrix_A;
01205 
01207         float64_t* observation_matrix_B;
01208 
01210         float64_t* transition_matrix_a;
01211 
01213         float64_t* initial_state_distribution_p;
01214 
01216         float64_t* end_state_distribution_q;        
01217 
01219         float64_t* observation_matrix_b;    
01220 
01222         int32_t iterations;
01223         int32_t iteration_count;
01224 
01226         float64_t epsilon;
01227         int32_t conv_it;
01228 
01230         float64_t all_pat_prob; 
01231 
01233         float64_t pat_prob; 
01234 
01236         float64_t mod_prob; 
01237 
01239         bool mod_prob_updated;  
01240 
01242         bool all_path_prob_updated; 
01243 
01245         int32_t path_deriv_dimension;
01246 
01248         bool path_deriv_updated;
01249 
01250         // true if model is using log likelihood
01251         bool loglikelihood;     
01252 
01253         // true->ok, false->error
01254         bool status;            
01255 
01256         // true->stolen from other HMMs, false->got own
01257         bool reused_caches;
01259 
01260 #ifdef USE_HMMPARALLEL_STRUCTURES
01261 
01262         float64_t** arrayN1 /*[parallel.get_num_threads()]*/ ;
01264         float64_t** arrayN2 /*[parallel.get_num_threads()]*/ ;
01265 #else //USE_HMMPARALLEL_STRUCTURES
01266 
01267         float64_t* arrayN1;
01269         float64_t* arrayN2;
01270 #endif //USE_HMMPARALLEL_STRUCTURES
01271 
01272 #ifdef USE_LOGSUMARRAY
01273 #ifdef USE_HMMPARALLEL_STRUCTURES
01274 
01275         float64_t** arrayS /*[parallel.get_num_threads()]*/;
01276 #else
01277 
01278         float64_t* arrayS;
01279 #endif // USE_HMMPARALLEL_STRUCTURES
01280 #endif // USE_LOGSUMARRAY
01281 
01282 #ifdef USE_HMMPARALLEL_STRUCTURES
01283 
01285         T_ALPHA_BETA* alpha_cache /*[parallel.get_num_threads()]*/ ;
01287         T_ALPHA_BETA* beta_cache /*[parallel.get_num_threads()]*/ ;
01288 
01290         T_STATES** states_per_observation_psi /*[parallel.get_num_threads()]*/ ;
01291 
01293         T_STATES** path /*[parallel.get_num_threads()]*/ ;
01294 
01296         bool* path_prob_updated /*[parallel.get_num_threads()]*/;
01297 
01299         int32_t* path_prob_dimension /*[parallel.get_num_threads()]*/ ; 
01300 
01301 #else //USE_HMMPARALLEL_STRUCTURES
01302 
01303         T_ALPHA_BETA alpha_cache;
01305         T_ALPHA_BETA beta_cache;
01306 
01308         T_STATES* states_per_observation_psi;
01309 
01311         T_STATES* path;
01312 
01314         bool path_prob_updated;
01315 
01317         int32_t path_prob_dimension;
01318 
01319 #endif //USE_HMMPARALLEL_STRUCTURES
01320 
01321 
01323         static const int32_t GOTN;
01325         static const int32_t GOTM;
01327         static const int32_t GOTO;
01329         static const int32_t GOTa;
01331         static const int32_t GOTb;
01333         static const int32_t GOTp;
01335         static const int32_t GOTq;
01336 
01338         static const int32_t GOTlearn_a;
01340         static const int32_t GOTlearn_b;
01342         static const int32_t GOTlearn_p;
01344         static const int32_t GOTlearn_q;
01346         static const int32_t GOTconst_a;
01348         static const int32_t GOTconst_b;
01350         static const int32_t GOTconst_p;
01352         static const int32_t GOTconst_q;
01353 
01354         public:
01359 
01361 inline float64_t state_probability(
01362     int32_t time, int32_t state, int32_t dimension)
01363 {
01364     return forward(time, state, dimension) + backward(time, state, dimension) - model_probability(dimension);
01365 }
01366 
01368 inline float64_t transition_probability(
01369     int32_t time, int32_t state_i, int32_t state_j, int32_t dimension)
01370 {
01371     return forward(time, state_i, dimension) + 
01372         backward(time+1, state_j, dimension) + 
01373         get_a(state_i,state_j) + get_b(state_j,p_observations->get_feature(dimension ,time+1)) - model_probability(dimension);
01374 }
01375 
01382 
01385 inline float64_t linear_model_derivative(
01386     T_STATES i, uint16_t j, int32_t dimension)
01387 {
01388     float64_t der=0;
01389 
01390     for (int32_t k=0; k<N; k++)
01391     {
01392         if (k!=i || p_observations->get_feature(dimension, k) != j)
01393             der+=get_b(k, p_observations->get_feature(dimension, k));
01394     }
01395 
01396     return der;
01397 }
01398 
01402 inline float64_t model_derivative_p(T_STATES i, int32_t dimension)
01403 {
01404     return backward(0,i,dimension)+get_b(i, p_observations->get_feature(dimension, 0));     
01405 }
01406 
01410 inline float64_t model_derivative_q(T_STATES i, int32_t dimension)
01411 {
01412     return forward(p_observations->get_vector_length(dimension)-1,i,dimension) ;
01413 }
01414 
01416 inline float64_t model_derivative_a(T_STATES i, T_STATES j, int32_t dimension)
01417 {
01418     float64_t sum=-CMath::INFTY;
01419     for (int32_t t=0; t<p_observations->get_vector_length(dimension)-1; t++)
01420         sum= CMath::logarithmic_sum(sum, forward(t, i, dimension) + backward(t+1, j, dimension) + get_b(j, p_observations->get_feature(dimension,t+1)));
01421 
01422     return sum;
01423 }
01424 
01425 
01427 inline float64_t model_derivative_b(T_STATES i, uint16_t j, int32_t dimension)
01428 {
01429     float64_t sum=-CMath::INFTY;
01430     for (int32_t t=0; t<p_observations->get_vector_length(dimension); t++)
01431     {
01432         if (p_observations->get_feature(dimension,t)==j)
01433             sum= CMath::logarithmic_sum(sum, forward(t,i,dimension)+backward(t,i,dimension)-get_b(i,p_observations->get_feature(dimension,t)));
01434     }
01435     //if (sum==-CMath::INFTY)
01436     // SG_DEBUG( "log derivative is -inf: dim=%i, state=%i, obs=%i\n",dimension, i, j) ;
01437     return sum;
01438 } 
01440 
01447 
01449 inline float64_t path_derivative_p(T_STATES i, int32_t dimension)
01450 {
01451     best_path(dimension);
01452     return (i==PATH(dimension)[0]) ? (exp(-get_p(PATH(dimension)[0]))) : (0) ;
01453 }
01454 
01456 inline float64_t path_derivative_q(T_STATES i, int32_t dimension)
01457 {
01458     best_path(dimension);
01459     return (i==PATH(dimension)[p_observations->get_vector_length(dimension)-1]) ? (exp(-get_q(PATH(dimension)[p_observations->get_vector_length(dimension)-1]))) : 0 ;
01460 }
01461 
01463 inline float64_t path_derivative_a(T_STATES i, T_STATES j, int32_t dimension)
01464 {
01465     prepare_path_derivative(dimension) ;
01466     return (get_A(i,j)==0) ? (0) : (get_A(i,j)*exp(-get_a(i,j))) ;
01467 }
01468 
01470 inline float64_t path_derivative_b(T_STATES i, uint16_t j, int32_t dimension)
01471 {
01472     prepare_path_derivative(dimension) ;
01473     return (get_B(i,j)==0) ? (0) : (get_B(i,j)*exp(-get_b(i,j))) ;
01474 } 
01475 
01477 
01478 
01479 protected:
01484 
01485     bool get_numbuffer(FILE* file, char* buffer, int32_t length);
01486 
01488     void open_bracket(FILE* file);
01489 
01491     void close_bracket(FILE* file);
01492 
01494     bool comma_or_space(FILE* file);
01495 
01497     inline void error(int32_t p_line, const char* str)
01498     {
01499         if (p_line)
01500             SG_ERROR( "error in line %d %s\n", p_line, str);
01501         else
01502             SG_ERROR( "error %s\n", str);
01503     }
01505 
01507     inline void prepare_path_derivative(int32_t dim)
01508     {
01509         if (path_deriv_updated && (path_deriv_dimension==dim))
01510             return ;
01511         int32_t i,j,t ;
01512         best_path(dim);
01513         //initialize with zeros
01514         for (i=0; i<N; i++)
01515         {
01516             for (j=0; j<N; j++)
01517                 set_A(i,j, 0);
01518             for (j=0; j<M; j++)
01519                 set_B(i,j, 0);
01520         }
01521 
01522         //counting occurences for A and B
01523         for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
01524         {
01525             set_A(PATH(dim)[t], PATH(dim)[t+1], get_A(PATH(dim)[t], PATH(dim)[t+1])+1);
01526             set_B(PATH(dim)[t], p_observations->get_feature(dim,t),  get_B(PATH(dim)[t], p_observations->get_feature(dim,t))+1);
01527         }
01528         set_B(PATH(dim)[p_observations->get_vector_length(dim)-1], p_observations->get_feature(dim,p_observations->get_vector_length(dim)-1),  get_B(PATH(dim)[p_observations->get_vector_length(dim)-1], p_observations->get_feature(dim,p_observations->get_vector_length(dim)-1)) + 1);
01529         path_deriv_dimension=dim ;
01530         path_deriv_updated=true ;
01531     } ;
01533 
01535     inline float64_t forward(int32_t time, int32_t state, int32_t dimension)
01536     {
01537         if (time<1)
01538             time=0;
01539 
01540         if (ALPHA_CACHE(dimension).table && (dimension==ALPHA_CACHE(dimension).dimension) && ALPHA_CACHE(dimension).updated)
01541         {
01542             if (time<p_observations->get_vector_length(dimension))
01543                 return ALPHA_CACHE(dimension).table[time*N+state];
01544             else
01545                 return ALPHA_CACHE(dimension).sum;
01546         }
01547         else
01548             return forward_comp(time, state, dimension) ;
01549     }
01550 
01552     inline float64_t backward(int32_t time, int32_t state, int32_t dimension)
01553     {
01554         if (BETA_CACHE(dimension).table && (dimension==BETA_CACHE(dimension).dimension) && (BETA_CACHE(dimension).updated))
01555         {
01556             if (time<0)
01557                 return BETA_CACHE(dimension).sum;
01558             if (time<p_observations->get_vector_length(dimension))
01559                 return BETA_CACHE(dimension).table[time*N+state];
01560             else
01561                 return -CMath::INFTY;
01562         }
01563         else
01564             return backward_comp(time, state, dimension) ;
01565     }
01566 
01567 };
01568 }
01569 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation