HMM.h

Go to the documentation of this file.

/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * Written (W) 1999-2009 Soeren Sonnenburg
 * Written (W) 1999-2008 Gunnar Raetsch
 * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
 */

#ifndef __CHMM_H__
#define __CHMM_H__

#include <shogun/mathematics/Math.h>
#include <shogun/lib/common.h>
#include <shogun/io/SGIO.h>
#include <shogun/lib/config.h>
#include <shogun/features/Features.h>
#include <shogun/features/StringFeatures.h>
#include <shogun/distributions/Distribution.h>

#include <stdio.h>

#ifdef USE_HMMPARALLEL
#define USE_HMMPARALLEL_STRUCTURES 1
#endif

namespace shogun
{
    class CFeatures;
    template <class ST> class CStringFeatures;

typedef float64_t T_ALPHA_BETA_TABLE;

#ifndef DOXYGEN_SHOULD_SKIP_THIS

struct T_ALPHA_BETA
{
    int32_t dimension;

    T_ALPHA_BETA_TABLE* table;

    bool updated;

    float64_t sum;
};
#endif // DOXYGEN_SHOULD_SKIP_THIS

#ifdef USE_BIGSTATES
typedef uint16_t T_STATES ;
#else
typedef uint8_t T_STATES ;
#endif
typedef T_STATES* P_STATES ;


enum BaumWelchViterbiType
{
    BW_NORMAL,
    BW_TRANS,
    BW_DEFINED,
    VIT_NORMAL,
    VIT_DEFINED
};


class Model
{
    public:
        Model();

        virtual ~Model();

        inline void sort_learn_a()
        {
            CMath::sort(learn_a,2) ;
        }

        inline void sort_learn_b()
        {
            CMath::sort(learn_b,2) ;
        }


        inline int32_t get_learn_a(int32_t line, int32_t column) const
        {
            return learn_a[line*2 + column];
        }

        inline int32_t get_learn_b(int32_t line, int32_t column) const
        {
            return learn_b[line*2 + column];
        }

        inline int32_t get_learn_p(int32_t offset) const
        {
            return learn_p[offset];
        }

        inline int32_t get_learn_q(int32_t offset) const
        {
            return learn_q[offset];
        }

        inline int32_t get_const_a(int32_t line, int32_t column) const
        {
            return const_a[line*2 + column];
        }

        inline int32_t get_const_b(int32_t line, int32_t column) const
        {
            return const_b[line*2 + column];
        }

        inline int32_t get_const_p(int32_t offset) const
        {
            return const_p[offset];
        }

        inline int32_t get_const_q(int32_t offset) const
        {
            return const_q[offset];
        }

        inline float64_t get_const_a_val(int32_t line) const
        {
            return const_a_val[line];
        }

        inline float64_t get_const_b_val(int32_t line) const
        {
            return const_b_val[line];
        }

        inline float64_t get_const_p_val(int32_t offset) const
        {
            return const_p_val[offset];
        }

        inline float64_t get_const_q_val(int32_t offset) const
        {
            return const_q_val[offset];
        }
#ifdef FIX_POS

        inline char get_fix_pos_state(int32_t pos, T_STATES state, T_STATES num_states)
        {
#ifdef HMM_DEBUG
            if ((pos<0)||(pos*num_states+state>65336))
                SG_DEBUG("index out of range in get_fix_pos_state(%i,%i,%i) \n", pos,state,num_states) ;
#endif
            return fix_pos_state[pos*num_states+state] ;
        }
#endif



        inline void set_learn_a(int32_t offset, int32_t value)
        {
            learn_a[offset]=value;
        }

        inline void set_learn_b(int32_t offset, int32_t value)
        {
            learn_b[offset]=value;
        }

        inline void set_learn_p(int32_t offset, int32_t value)
        {
            learn_p[offset]=value;
        }

        inline void set_learn_q(int32_t offset, int32_t value)
        {
            learn_q[offset]=value;
        }

        inline void set_const_a(int32_t offset, int32_t value)
        {
            const_a[offset]=value;
        }

        inline void set_const_b(int32_t offset, int32_t value)
        {
            const_b[offset]=value;
        }

        inline void set_const_p(int32_t offset, int32_t value)
        {
            const_p[offset]=value;
        }

        inline void set_const_q(int32_t offset, int32_t value)
        {
            const_q[offset]=value;
        }

        inline void set_const_a_val(int32_t offset, float64_t value)
        {
            const_a_val[offset]=value;
        }

        inline void set_const_b_val(int32_t offset, float64_t value)
        {
            const_b_val[offset]=value;
        }

        inline void set_const_p_val(int32_t offset, float64_t value)
        {
            const_p_val[offset]=value;
        }

        inline void set_const_q_val(int32_t offset, float64_t value)
        {
            const_q_val[offset]=value;
        }
#ifdef FIX_POS

        inline void set_fix_pos_state(
            int32_t pos, T_STATES state, T_STATES num_states, char value)
        {
#ifdef HMM_DEBUG
            if ((pos<0)||(pos*num_states+state>65336))
                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) ;
#endif
            fix_pos_state[pos*num_states+state]=value;
            if (value==FIX_ALLOWED)
                for (int32_t i=0; i<num_states; i++)
                    if (get_fix_pos_state(pos,i,num_states)==FIX_DEFAULT)
                        set_fix_pos_state(pos,i,num_states,FIX_DISALLOWED) ;
        }

        const static char FIX_DISALLOWED ;

        const static char FIX_ALLOWED ;

        const static char FIX_DEFAULT ;

        const static float64_t DISALLOWED_PENALTY ;
#endif
    protected:

        int32_t* learn_a;

        int32_t* learn_b;

        int32_t* learn_p;

        int32_t* learn_q;


        int32_t* const_a;

        int32_t* const_b;

        int32_t* const_p;

        int32_t* const_q;


        float64_t* const_a_val;

        float64_t* const_b_val;

        float64_t* const_p_val;

        float64_t* const_q_val;

#ifdef FIX_POS

        char* fix_pos_state;
#endif

};


class CHMM : public CDistribution
{
    private:

        T_STATES trans_list_len ;
        T_STATES **trans_list_forward  ;
        T_STATES *trans_list_forward_cnt  ;
        float64_t **trans_list_forward_val ;
        T_STATES **trans_list_backward  ;
        T_STATES *trans_list_backward_cnt  ;
        bool mem_initialized ;

#ifdef USE_HMMPARALLEL_STRUCTURES

        struct S_DIM_THREAD_PARAM
        {
            CHMM* hmm;
            int32_t dim;
            float64_t prob_sum;
        };

        struct S_BW_THREAD_PARAM
        {
            CHMM* hmm;
            int32_t dim_start;
            int32_t dim_stop;

            float64_t ret;

            float64_t* p_buf;
            float64_t* q_buf;
            float64_t* a_buf;
            float64_t* b_buf;
        };

        inline T_ALPHA_BETA & ALPHA_CACHE(int32_t dim) {
            return alpha_cache[dim%parallel->get_num_threads()] ; } ;
        inline T_ALPHA_BETA & BETA_CACHE(int32_t dim) {
            return beta_cache[dim%parallel->get_num_threads()] ; } ;
#ifdef USE_LOGSUMARRAY
        inline float64_t* ARRAYS(int32_t dim) {
            return arrayS[dim%parallel->get_num_threads()] ; } ;
#endif
        inline float64_t* ARRAYN1(int32_t dim) {
            return arrayN1[dim%parallel->get_num_threads()] ; } ;
        inline float64_t* ARRAYN2(int32_t dim) {
            return arrayN2[dim%parallel->get_num_threads()] ; } ;
        inline T_STATES* STATES_PER_OBSERVATION_PSI(int32_t dim) {
            return states_per_observation_psi[dim%parallel->get_num_threads()] ; } ;
        inline const T_STATES* STATES_PER_OBSERVATION_PSI(int32_t dim) const {
            return states_per_observation_psi[dim%parallel->get_num_threads()] ; } ;
        inline T_STATES* PATH(int32_t dim) {
            return path[dim%parallel->get_num_threads()] ; } ;
        inline bool & PATH_PROB_UPDATED(int32_t dim) {
            return path_prob_updated[dim%parallel->get_num_threads()] ; } ;
        inline int32_t & PATH_PROB_DIMENSION(int32_t dim) {
            return path_prob_dimension[dim%parallel->get_num_threads()] ; } ;
#else
        inline T_ALPHA_BETA & ALPHA_CACHE(int32_t /*dim*/) {
            return alpha_cache ; } ;
        inline T_ALPHA_BETA & BETA_CACHE(int32_t /*dim*/) {
            return beta_cache ; } ;
#ifdef USE_LOGSUMARRAY
        inline float64_t* ARRAYS(int32_t dim) {
            return arrayS ; } ;
#endif
        inline float64_t* ARRAYN1(int32_t /*dim*/) {
            return arrayN1 ; } ;
        inline float64_t* ARRAYN2(int32_t /*dim*/) {
            return arrayN2 ; } ;
        inline T_STATES* STATES_PER_OBSERVATION_PSI(int32_t /*dim*/) {
            return states_per_observation_psi ; } ;
        inline const T_STATES* STATES_PER_OBSERVATION_PSI(int32_t /*dim*/) const {
            return states_per_observation_psi ; } ;
        inline T_STATES* PATH(int32_t /*dim*/) {
            return path ; } ;
        inline bool & PATH_PROB_UPDATED(int32_t /*dim*/) {
            return path_prob_updated ; } ;
        inline int32_t & PATH_PROB_DIMENSION(int32_t /*dim*/) {
            return path_prob_dimension ; } ;
#endif

        bool converged(float64_t x, float64_t y);

    public:
        CHMM();

        CHMM(
            int32_t N, int32_t M, Model* model, float64_t PSEUDO);
        CHMM(
            CStringFeatures<uint16_t>* obs, int32_t N, int32_t M,
            float64_t PSEUDO);
        CHMM(
            int32_t N, float64_t* p, float64_t* q, float64_t* a);
        CHMM(
            int32_t N, float64_t* p, float64_t* q, int32_t num_trans,
            float64_t* a_trans);

        CHMM(FILE* model_file, float64_t PSEUDO);

        CHMM(CHMM* h);

        virtual ~CHMM();

        virtual bool train(CFeatures* data=NULL);
        virtual int32_t get_num_model_parameters() { return N*(N+M+2); }
        virtual float64_t get_log_model_parameter(int32_t num_param);
        virtual float64_t get_log_derivative(int32_t num_param, int32_t num_example);
        virtual float64_t get_log_likelihood_example(int32_t num_example)
        {
            return model_probability(num_example);
        }

        bool initialize(Model* model, float64_t PSEUDO, FILE* model_file=NULL);

        bool alloc_state_dependend_arrays();

        void free_state_dependend_arrays();

        float64_t forward_comp(int32_t time, int32_t state, int32_t dimension);
        float64_t forward_comp_old(
            int32_t time, int32_t state, int32_t dimension);

        float64_t backward_comp(int32_t time, int32_t state, int32_t dimension);
        float64_t backward_comp_old(
            int32_t time, int32_t state, int32_t dimension);

        float64_t best_path(int32_t dimension);
        inline uint16_t get_best_path_state(int32_t dim, int32_t t)
        {
            ASSERT(PATH(dim));
            return PATH(dim)[t];
        }

        float64_t model_probability_comp() ;

        inline float64_t model_probability(int32_t dimension=-1)
        {
            //for faster calculation cache model probability
            if (dimension==-1)
            {
                if (mod_prob_updated)
                    return mod_prob/p_observations->get_num_vectors();
                else
                    return model_probability_comp()/p_observations->get_num_vectors();
            }
            else
                return forward(p_observations->get_vector_length(dimension), 0, dimension);
        }

        inline float64_t linear_model_probability(int32_t dimension)
        {
            float64_t lik=0;
            int32_t len=0;
            bool free_vec;
            uint16_t* o=p_observations->get_feature_vector(dimension, len, free_vec);
            float64_t* obs_b=observation_matrix_b;

            ASSERT(N==len);

            for (int32_t i=0; i<N; i++)
            {
                lik+=obs_b[*o++];
                obs_b+=M;
            }
            p_observations->free_feature_vector(o, dimension, free_vec);
            return lik;

            // sorry, the above code is the speed optimized version of :
            /*  float64_t lik=0;

                for (int32_t i=0; i<N; i++)
                lik+=get_b(i, p_observations->get_feature(dimension, i));
                return lik;
                */
            // : that
        }


        inline bool set_iterations(int32_t num) { iterations=num; return true; }
        inline int32_t get_iterations() { return iterations; }
        inline bool set_epsilon (float64_t eps) { epsilon=eps; return true; }
        inline float64_t get_epsilon() { return epsilon; }

        bool baum_welch_viterbi_train(BaumWelchViterbiType type);

        void estimate_model_baum_welch(CHMM* train);
        void estimate_model_baum_welch_trans(CHMM* train);

#ifdef USE_HMMPARALLEL_STRUCTURES
        void ab_buf_comp(
            float64_t* p_buf, float64_t* q_buf, float64_t* a_buf,
            float64_t* b_buf, int32_t dim) ;
#else
        void estimate_model_baum_welch_old(CHMM* train);
#endif

        void estimate_model_baum_welch_defined(CHMM* train);

        void estimate_model_viterbi(CHMM* train);

        void estimate_model_viterbi_defined(CHMM* train);


        bool linear_train(bool right_align=false);

        bool permutation_entropy(int32_t window_width, int32_t sequence_number);

        void output_model(bool verbose=false);

        void output_model_defined(bool verbose=false);



        void normalize(bool keep_dead_states=false);

        void add_states(int32_t num_states, float64_t default_val=0);

        bool append_model(
            CHMM* append_model, float64_t* cur_out, float64_t* app_out);

        bool append_model(CHMM* append_model);

        void chop(float64_t value);

        void convert_to_log();

        void init_model_random();

        void init_model_defined();

        void clear_model();

        void clear_model_defined();

        void copy_model(CHMM* l);

        void invalidate_model();

        inline bool get_status() const
        {
            return status;
        }

        inline float64_t get_pseudo() const
        {
            return PSEUDO ;
        }

        inline void set_pseudo(float64_t pseudo)
        {
            PSEUDO=pseudo ;
        }

#ifdef USE_HMMPARALLEL_STRUCTURES
        static void* bw_dim_prefetch(void * params);
        static void* bw_single_dim_prefetch(void * params);
        static void* vit_dim_prefetch(void * params);
#endif

#ifdef FIX_POS

        inline bool set_fix_pos_state(int32_t pos, T_STATES state, char value)
        {
            if (!model)
                return false ;
            model->set_fix_pos_state(pos, state, N, value) ;
            return true ;
        } ;
#endif


        void set_observations(CStringFeatures<uint16_t>* obs, CHMM* hmm=NULL);

        void set_observation_nocache(CStringFeatures<uint16_t>* obs);

        inline CStringFeatures<uint16_t>* get_observations()
        {
            SG_REF(p_observations);
            return p_observations;
        }

        bool load_definitions(FILE* file, bool verbose, bool initialize=true);

        bool load_model(FILE* file);

        bool save_model(FILE* file);

        bool save_model_derivatives(FILE* file);

        bool save_model_derivatives_bin(FILE* file);

        bool save_model_bin(FILE* file);

        bool check_model_derivatives() ;
        bool check_model_derivatives_combined() ;

        T_STATES* get_path(int32_t dim, float64_t& prob);

        bool save_path(FILE* file);

        bool save_path_derivatives(FILE* file);

        bool save_path_derivatives_bin(FILE* file);

#ifdef USE_HMMDEBUG

        bool check_path_derivatives() ;
#endif //USE_HMMDEBUG

        bool save_likelihood_bin(FILE* file);

        bool save_likelihood(FILE* file);


        inline T_STATES get_N() const { return N ; }

        inline int32_t get_M() const { return M ; }

        inline void set_q(T_STATES offset, float64_t value)
        {
#ifdef HMM_DEBUG
            if (offset>=N)
                SG_DEBUG("index out of range in set_q(%i,%e) [%i]\n", offset,value,N) ;
#endif
            end_state_distribution_q[offset]=value;
        }

        inline void set_p(T_STATES offset, float64_t value)
        {
#ifdef HMM_DEBUG
            if (offset>=N)
                SG_DEBUG("index out of range in set_p(%i,.) [%i]\n", offset,N) ;
#endif
            initial_state_distribution_p[offset]=value;
        }

        inline void set_A(T_STATES line_, T_STATES column, float64_t value)
        {
#ifdef HMM_DEBUG
            if ((line_>N)||(column>N))
                SG_DEBUG("index out of range in set_A(%i,%i,.) [%i,%i]\n",line_,column,N,N) ;
#endif
            transition_matrix_A[line_+column*N]=value;
        }

        inline void set_a(T_STATES line_, T_STATES column, float64_t value)
        {
#ifdef HMM_DEBUG
            if ((line_>N)||(column>N))
                SG_DEBUG("index out of range in set_a(%i,%i,.) [%i,%i]\n",line_,column,N,N) ;
#endif
            transition_matrix_a[line_+column*N]=value; // look also best_path!
        }

        inline void set_B(T_STATES line_, uint16_t column, float64_t value)
        {
#ifdef HMM_DEBUG
            if ((line_>=N)||(column>=M))
                SG_DEBUG("index out of range in set_B(%i,%i) [%i,%i]\n", line_, column,N,M) ;
#endif
            observation_matrix_B[line_*M+column]=value;
        }

        inline void set_b(T_STATES line_, uint16_t column, float64_t value)
        {
#ifdef HMM_DEBUG
            if ((line_>=N)||(column>=M))
                SG_DEBUG("index out of range in set_b(%i,%i) [%i,%i]\n", line_, column,N,M) ;
#endif
            observation_matrix_b[line_*M+column]=value;
        }

        inline void set_psi(
            int32_t time, T_STATES state, T_STATES value, int32_t dimension)
        {
#ifdef HMM_DEBUG
            if ((time>=p_observations->get_max_vector_length())||(state>N))
                SG_DEBUG("index out of range in set_psi(%i,%i,.) [%i,%i]\n",time,state,p_observations->get_max_vector_length(),N) ;
#endif
            STATES_PER_OBSERVATION_PSI(dimension)[time*N+state]=value;
        }

        inline float64_t get_q(T_STATES offset) const
        {
#ifdef HMM_DEBUG
            if (offset>=N)
                SG_DEBUG("index out of range in %e=get_q(%i) [%i]\n", end_state_distribution_q[offset],offset,N) ;
#endif
            return end_state_distribution_q[offset];
        }

        inline float64_t get_p(T_STATES offset) const
        {
#ifdef HMM_DEBUG
            if (offset>=N)
                SG_DEBUG("index out of range in get_p(%i,.) [%i]\n", offset,N) ;
#endif
            return initial_state_distribution_p[offset];
        }

        inline float64_t get_A(T_STATES line_, T_STATES column) const
        {
#ifdef HMM_DEBUG
            if ((line_>N)||(column>N))
                SG_DEBUG("index out of range in get_A(%i,%i) [%i,%i]\n",line_,column,N,N) ;
#endif
            return transition_matrix_A[line_+column*N];
        }

        inline float64_t get_a(T_STATES line_, T_STATES column) const
        {
#ifdef HMM_DEBUG
            if ((line_>N)||(column>N))
                SG_DEBUG("index out of range in get_a(%i,%i) [%i,%i]\n",line_,column,N,N) ;
#endif
            return transition_matrix_a[line_+column*N]; // look also best_path()!
        }

        inline float64_t get_B(T_STATES line_, uint16_t column) const
        {
#ifdef HMM_DEBUG
            if ((line_>=N)||(column>=M))
                SG_DEBUG("index out of range in get_B(%i,%i) [%i,%i]\n", line_, column,N,M) ;
#endif
            return observation_matrix_B[line_*M+column];
        }

        inline float64_t get_b(T_STATES line_, uint16_t column) const
        {
#ifdef HMM_DEBUG
            if ((line_>=N)||(column>=M))
                SG_DEBUG("index out of range in get_b(%i,%i) [%i,%i]\n", line_, column,N,M) ;
#endif
            //SG_PRINT("idx %d\n", line_*M+column);
            return observation_matrix_b[line_*M+column];
        }

        inline T_STATES get_psi(
            int32_t time, T_STATES state, int32_t dimension) const
        {
#ifdef HMM_DEBUG
            if ((time>=p_observations->get_max_vector_length())||(state>N))
                SG_DEBUG("index out of range in get_psi(%i,%i) [%i,%i]\n",time,state,p_observations->get_max_vector_length(),N) ;
#endif
            return STATES_PER_OBSERVATION_PSI(dimension)[time*N+state];
        }


        virtual const char* get_name() const { return "HMM"; }

    protected:

        int32_t M;

        int32_t N;

        float64_t PSEUDO;

        // line number during processing input files
        int32_t line;

        CStringFeatures<uint16_t>* p_observations;

        //train definition for HMM
        Model* model;

        float64_t* transition_matrix_A;

        float64_t* observation_matrix_B;

        float64_t* transition_matrix_a;

        float64_t* initial_state_distribution_p;

        float64_t* end_state_distribution_q;

        float64_t* observation_matrix_b;

        int32_t iterations;
        int32_t iteration_count;

        float64_t epsilon;
        int32_t conv_it;

        float64_t all_pat_prob;

        float64_t pat_prob;

        float64_t mod_prob;

        bool mod_prob_updated;

        bool all_path_prob_updated;

        int32_t path_deriv_dimension;

        bool path_deriv_updated;

        // true if model is using log likelihood
        bool loglikelihood;

        // true->ok, false->error
        bool status;

        // true->stolen from other HMMs, false->got own
        bool reused_caches;

#ifdef USE_HMMPARALLEL_STRUCTURES

        float64_t** arrayN1 /*[parallel.get_num_threads()]*/ ;
        float64_t** arrayN2 /*[parallel.get_num_threads()]*/ ;
#else //USE_HMMPARALLEL_STRUCTURES

        float64_t* arrayN1;
        float64_t* arrayN2;
#endif //USE_HMMPARALLEL_STRUCTURES

#ifdef USE_LOGSUMARRAY
#ifdef USE_HMMPARALLEL_STRUCTURES

        float64_t** arrayS /*[parallel.get_num_threads()]*/;
#else

        float64_t* arrayS;
#endif // USE_HMMPARALLEL_STRUCTURES
#endif // USE_LOGSUMARRAY

#ifdef USE_HMMPARALLEL_STRUCTURES

        T_ALPHA_BETA* alpha_cache /*[parallel.get_num_threads()]*/ ;
        T_ALPHA_BETA* beta_cache /*[parallel.get_num_threads()]*/ ;

        T_STATES** states_per_observation_psi /*[parallel.get_num_threads()]*/ ;

        T_STATES** path /*[parallel.get_num_threads()]*/ ;

        bool* path_prob_updated /*[parallel.get_num_threads()]*/;

        int32_t* path_prob_dimension /*[parallel.get_num_threads()]*/ ;

#else //USE_HMMPARALLEL_STRUCTURES

        T_ALPHA_BETA alpha_cache;
        T_ALPHA_BETA beta_cache;

        T_STATES* states_per_observation_psi;

        T_STATES* path;

        bool path_prob_updated;

        int32_t path_prob_dimension;

#endif //USE_HMMPARALLEL_STRUCTURES


        static const int32_t GOTN;
        static const int32_t GOTM;
        static const int32_t GOTO;
        static const int32_t GOTa;
        static const int32_t GOTb;
        static const int32_t GOTp;
        static const int32_t GOTq;

        static const int32_t GOTlearn_a;
        static const int32_t GOTlearn_b;
        static const int32_t GOTlearn_p;
        static const int32_t GOTlearn_q;
        static const int32_t GOTconst_a;
        static const int32_t GOTconst_b;
        static const int32_t GOTconst_p;
        static const int32_t GOTconst_q;

        public:

inline float64_t state_probability(
    int32_t time, int32_t state, int32_t dimension)
{
    return forward(time, state, dimension) + backward(time, state, dimension) - model_probability(dimension);
}

inline float64_t transition_probability(
    int32_t time, int32_t state_i, int32_t state_j, int32_t dimension)
{
    return forward(time, state_i, dimension) +
        backward(time+1, state_j, dimension) +
        get_a(state_i,state_j) + get_b(state_j,p_observations->get_feature(dimension ,time+1)) - model_probability(dimension);
}


inline float64_t linear_model_derivative(
    T_STATES i, uint16_t j, int32_t dimension)
{
    float64_t der=0;

    for (int32_t k=0; k<N; k++)
    {
        if (k!=i || p_observations->get_feature(dimension, k) != j)
            der+=get_b(k, p_observations->get_feature(dimension, k));
    }

    return der;
}

inline float64_t model_derivative_p(T_STATES i, int32_t dimension)
{
    return backward(0,i,dimension)+get_b(i, p_observations->get_feature(dimension, 0));
}

inline float64_t model_derivative_q(T_STATES i, int32_t dimension)
{
    return forward(p_observations->get_vector_length(dimension)-1,i,dimension) ;
}

inline float64_t model_derivative_a(T_STATES i, T_STATES j, int32_t dimension)
{
    float64_t sum=-CMath::INFTY;
    for (int32_t t=0; t<p_observations->get_vector_length(dimension)-1; t++)
        sum= CMath::logarithmic_sum(sum, forward(t, i, dimension) + backward(t+1, j, dimension) + get_b(j, p_observations->get_feature(dimension,t+1)));

    return sum;
}


inline float64_t model_derivative_b(T_STATES i, uint16_t j, int32_t dimension)
{
    float64_t sum=-CMath::INFTY;
    for (int32_t t=0; t<p_observations->get_vector_length(dimension); t++)
    {
        if (p_observations->get_feature(dimension,t)==j)
            sum= CMath::logarithmic_sum(sum, forward(t,i,dimension)+backward(t,i,dimension)-get_b(i,p_observations->get_feature(dimension,t)));
    }
    //if (sum==-CMath::INFTY)
    // SG_DEBUG( "log derivative is -inf: dim=%i, state=%i, obs=%i\n",dimension, i, j) ;
    return sum;
}


inline float64_t path_derivative_p(T_STATES i, int32_t dimension)
{
    best_path(dimension);
    return (i==PATH(dimension)[0]) ? (exp(-get_p(PATH(dimension)[0]))) : (0) ;
}

inline float64_t path_derivative_q(T_STATES i, int32_t dimension)
{
    best_path(dimension);
    return (i==PATH(dimension)[p_observations->get_vector_length(dimension)-1]) ? (exp(-get_q(PATH(dimension)[p_observations->get_vector_length(dimension)-1]))) : 0 ;
}

inline float64_t path_derivative_a(T_STATES i, T_STATES j, int32_t dimension)
{
    prepare_path_derivative(dimension) ;
    return (get_A(i,j)==0) ? (0) : (get_A(i,j)*exp(-get_a(i,j))) ;
}

inline float64_t path_derivative_b(T_STATES i, uint16_t j, int32_t dimension)
{
    prepare_path_derivative(dimension) ;
    return (get_B(i,j)==0) ? (0) : (get_B(i,j)*exp(-get_b(i,j))) ;
}



protected:

    bool get_numbuffer(FILE* file, char* buffer, int32_t length);

    void open_bracket(FILE* file);

    void close_bracket(FILE* file);

    bool comma_or_space(FILE* file);

    inline void error(int32_t p_line, const char* str)
    {
        if (p_line)
            SG_ERROR( "error in line %d %s\n", p_line, str);
        else
            SG_ERROR( "error %s\n", str);
    }

    inline void prepare_path_derivative(int32_t dim)
    {
        if (path_deriv_updated && (path_deriv_dimension==dim))
            return ;
        int32_t i,j,t ;
        best_path(dim);
        //initialize with zeros
        for (i=0; i<N; i++)
        {
            for (j=0; j<N; j++)
                set_A(i,j, 0);
            for (j=0; j<M; j++)
                set_B(i,j, 0);
        }

        //counting occurences for A and B
        for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
        {
            set_A(PATH(dim)[t], PATH(dim)[t+1], get_A(PATH(dim)[t], PATH(dim)[t+1])+1);
            set_B(PATH(dim)[t], p_observations->get_feature(dim,t),  get_B(PATH(dim)[t], p_observations->get_feature(dim,t))+1);
        }
        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);
        path_deriv_dimension=dim ;
        path_deriv_updated=true ;
    } ;

    inline float64_t forward(int32_t time, int32_t state, int32_t dimension)
    {
        if (time<1)
            time=0;

        if (ALPHA_CACHE(dimension).table && (dimension==ALPHA_CACHE(dimension).dimension) && ALPHA_CACHE(dimension).updated)
        {
            if (time<p_observations->get_vector_length(dimension))
                return ALPHA_CACHE(dimension).table[time*N+state];
            else
                return ALPHA_CACHE(dimension).sum;
        }
        else
            return forward_comp(time, state, dimension) ;
    }

    inline float64_t backward(int32_t time, int32_t state, int32_t dimension)
    {
        if (BETA_CACHE(dimension).table && (dimension==BETA_CACHE(dimension).dimension) && (BETA_CACHE(dimension).updated))
        {
            if (time<0)
                return BETA_CACHE(dimension).sum;
            if (time<p_observations->get_vector_length(dimension))
                return BETA_CACHE(dimension).table[time*N+state];
            else
                return -CMath::INFTY;
        }
        else
            return backward_comp(time, state, dimension) ;
    }

};
}
#endif