HMM.cpp

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
 */
#include <shogun/distributions/HMM.h>
#include <shogun/mathematics/Math.h>
#include <shogun/io/SGIO.h>
#include <shogun/lib/config.h>
#include <shogun/lib/Signal.h>
#include <shogun/base/Parallel.h>
#include <shogun/features/StringFeatures.h>
#include <shogun/features/Alphabet.h>

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <ctype.h>

#define VAL_MACRO log((default_value == 0) ? (CMath::random(MIN_RAND, MAX_RAND)) : default_value)
#define ARRAY_SIZE 65336

using namespace shogun;

// Construction/Destruction

const int32_t CHMM::GOTN= (1<<1);
const int32_t CHMM::GOTM= (1<<2);
const int32_t CHMM::GOTO= (1<<3);
const int32_t CHMM::GOTa= (1<<4);
const int32_t CHMM::GOTb= (1<<5);
const int32_t CHMM::GOTp= (1<<6);
const int32_t CHMM::GOTq= (1<<7);

const int32_t CHMM::GOTlearn_a= (1<<1);
const int32_t CHMM::GOTlearn_b= (1<<2);
const int32_t CHMM::GOTlearn_p= (1<<3);
const int32_t CHMM::GOTlearn_q= (1<<4);
const int32_t CHMM::GOTconst_a= (1<<5);
const int32_t CHMM::GOTconst_b= (1<<6);
const int32_t CHMM::GOTconst_p= (1<<7);
const int32_t CHMM::GOTconst_q= (1<<8);

enum E_STATE
{
    INITIAL,
    ARRAYs,
    GET_N,
    GET_M,
    GET_a,
    GET_b,
    GET_p,
    GET_q,
    GET_learn_a,
    GET_learn_b,
    GET_learn_p,
    GET_learn_q,
    GET_const_a,
    GET_const_b,
    GET_const_p,
    GET_const_q,
    COMMENT,
    END
};


#ifdef FIX_POS
const char Model::FIX_DISALLOWED=0 ;
const char Model::FIX_ALLOWED=1 ;
const char Model::FIX_DEFAULT=-1 ;
const float64_t Model::DISALLOWED_PENALTY=CMath::ALMOST_NEG_INFTY ;
#endif

Model::Model()
{
    const_a=SG_MALLOC(int, ARRAY_SIZE);             
    const_b=SG_MALLOC(int, ARRAY_SIZE);
    const_p=SG_MALLOC(int, ARRAY_SIZE);
    const_q=SG_MALLOC(int, ARRAY_SIZE);
    const_a_val=SG_MALLOC(float64_t, ARRAY_SIZE);           
    const_b_val=SG_MALLOC(float64_t, ARRAY_SIZE);
    const_p_val=SG_MALLOC(float64_t, ARRAY_SIZE);
    const_q_val=SG_MALLOC(float64_t, ARRAY_SIZE);


    learn_a=SG_MALLOC(int, ARRAY_SIZE);
    learn_b=SG_MALLOC(int, ARRAY_SIZE);
    learn_p=SG_MALLOC(int, ARRAY_SIZE);
    learn_q=SG_MALLOC(int, ARRAY_SIZE);

#ifdef FIX_POS
    fix_pos_state = SG_MALLOC(char, ARRAY_SIZE);
#endif
    for (int32_t i=0; i<ARRAY_SIZE; i++)
    {
        const_a[i]=-1 ;
        const_b[i]=-1 ;
        const_p[i]=-1 ;
        const_q[i]=-1 ;
        const_a_val[i]=1.0 ;
        const_b_val[i]=1.0 ;
        const_p_val[i]=1.0 ;
        const_q_val[i]=1.0 ;
        learn_a[i]=-1 ;
        learn_b[i]=-1 ;
        learn_p[i]=-1 ;
        learn_q[i]=-1 ;
#ifdef FIX_POS
        fix_pos_state[i] = FIX_DEFAULT ;
#endif
    } ;
}

Model::~Model()
{
    SG_FREE(const_a);
    SG_FREE(const_b);
    SG_FREE(const_p);
    SG_FREE(const_q);
    SG_FREE(const_a_val);
    SG_FREE(const_b_val);
    SG_FREE(const_p_val);
    SG_FREE(const_q_val);

    SG_FREE(learn_a);
    SG_FREE(learn_b);
    SG_FREE(learn_p);
    SG_FREE(learn_q);

#ifdef FIX_POS
    SG_FREE(fix_pos_state);
#endif

}

CHMM::CHMM()
{
    N=0;
    M=0;
    model=NULL;
    status=false;
    p_observations=NULL;
    trans_list_forward_cnt=NULL;
    trans_list_backward_cnt=NULL;
    trans_list_forward=NULL;
    trans_list_backward=NULL;
    trans_list_forward_val=NULL;
    iterations=150;
    epsilon=1e-4;
    conv_it=5;
    path=NULL;
    arrayN1=NULL;
    arrayN2=NULL;
    reused_caches=false;
    transition_matrix_a=NULL;
    observation_matrix_b=NULL;
    initial_state_distribution_p=NULL;
    end_state_distribution_q=NULL;
#ifdef USE_LOGSUMARRAY
    arrayS = NULL;
#endif
#ifdef USE_HMMPARALLEL_STRUCTURES
    this->alpha_cache=NULL;
    this->beta_cache=NULL;
    path_prob_updated = NULL;
    path_prob_dimension = NULL;
#else
    this->alpha_cache.table=NULL;
    this->beta_cache.table=NULL;
    this->alpha_cache.dimension=0;
    this->beta_cache.dimension=0;
#endif
    states_per_observation_psi=NULL;
    mem_initialized = false;
}

CHMM::CHMM(CHMM* h)
: CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
{
#ifdef USE_HMMPARALLEL_STRUCTURES
    SG_INFO("hmm is using %i separate tables\n",  parallel->get_num_threads())
#endif

    this->N=h->get_N();
    this->M=h->get_M();
    status=initialize_hmm(NULL, h->get_pseudo());
    this->copy_model(h);
    set_observations(h->p_observations);
}

CHMM::CHMM(int32_t p_N, int32_t p_M, Model* p_model, float64_t p_PSEUDO)
: CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
{
    this->N=p_N;
    this->M=p_M;
    model=NULL ;

#ifdef USE_HMMPARALLEL_STRUCTURES
    SG_INFO("hmm is using %i separate tables\n",  parallel->get_num_threads())
#endif

    status=initialize_hmm(p_model, p_PSEUDO);
}

CHMM::CHMM(
    CStringFeatures<uint16_t>* obs, int32_t p_N, int32_t p_M,
    float64_t p_PSEUDO)
: CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
{
    this->N=p_N;
    this->M=p_M;
    model=NULL ;

#ifdef USE_HMMPARALLEL_STRUCTURES
    SG_INFO("hmm is using %i separate tables\n",  parallel->get_num_threads())
#endif

    initialize_hmm(model, p_PSEUDO);
    set_observations(obs);
}

CHMM::CHMM(int32_t p_N, float64_t* p, float64_t* q, float64_t* a)
: CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
{
    this->N=p_N;
    this->M=0;
    model=NULL ;

    trans_list_forward = NULL ;
    trans_list_forward_cnt = NULL ;
    trans_list_forward_val = NULL ;
    trans_list_backward = NULL ;
    trans_list_backward_cnt = NULL ;
    trans_list_len = 0 ;
    mem_initialized = false ;

    this->transition_matrix_a=NULL;
    this->observation_matrix_b=NULL;
    this->initial_state_distribution_p=NULL;
    this->end_state_distribution_q=NULL;
    this->p_observations=NULL;
    this->reused_caches=false;

#ifdef USE_HMMPARALLEL_STRUCTURES
    this->alpha_cache=NULL;
    this->beta_cache=NULL;
#else
    this->alpha_cache.table=NULL;
    this->beta_cache.table=NULL;
    this->alpha_cache.dimension=0;
    this->beta_cache.dimension=0;
#endif

    this->states_per_observation_psi=NULL ;
    this->path=NULL;
    arrayN1=NULL ;
    arrayN2=NULL ;

    this->loglikelihood=false;
    mem_initialized = true ;

    transition_matrix_a=a ;
    observation_matrix_b=NULL ;
    initial_state_distribution_p=p ;
    end_state_distribution_q=q ;
    transition_matrix_A=NULL ;
    observation_matrix_B=NULL ;

//  this->invalidate_model();
}

CHMM::CHMM(
    int32_t p_N, float64_t* p, float64_t* q, int32_t num_trans,
    float64_t* a_trans)
: CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
{
    model=NULL ;

    this->N=p_N;
    this->M=0;

    trans_list_forward = NULL ;
    trans_list_forward_cnt = NULL ;
    trans_list_forward_val = NULL ;
    trans_list_backward = NULL ;
    trans_list_backward_cnt = NULL ;
    trans_list_len = 0 ;
    mem_initialized = false ;

    this->transition_matrix_a=NULL;
    this->observation_matrix_b=NULL;
    this->initial_state_distribution_p=NULL;
    this->end_state_distribution_q=NULL;
    this->p_observations=NULL;
    this->reused_caches=false;

#ifdef USE_HMMPARALLEL_STRUCTURES
    this->alpha_cache=NULL;
    this->beta_cache=NULL;
#else
    this->alpha_cache.table=NULL;
    this->beta_cache.table=NULL;
    this->alpha_cache.dimension=0;
    this->beta_cache.dimension=0;
#endif

    this->states_per_observation_psi=NULL ;
    this->path=NULL;
    arrayN1=NULL ;
    arrayN2=NULL ;

    this->loglikelihood=false;
    mem_initialized = true ;

    trans_list_forward_cnt=NULL ;
    trans_list_len = N ;
    trans_list_forward = SG_MALLOC(T_STATES*, N);
    trans_list_forward_val = SG_MALLOC(float64_t*, N);
    trans_list_forward_cnt = SG_MALLOC(T_STATES, N);

    int32_t start_idx=0;
    for (int32_t j=0; j<N; j++)
    {
        int32_t old_start_idx=start_idx;

        while (start_idx<num_trans && a_trans[start_idx+num_trans]==j)
        {
            start_idx++;

            if (start_idx>1 && start_idx<num_trans)
                ASSERT(a_trans[start_idx+num_trans-1]<=
                    a_trans[start_idx+num_trans]);
        }

        if (start_idx>1 && start_idx<num_trans)
            ASSERT(a_trans[start_idx+num_trans-1]<=
                a_trans[start_idx+num_trans]);

        int32_t len=start_idx-old_start_idx;
        ASSERT(len>=0)

        trans_list_forward_cnt[j] = 0 ;

        if (len>0)
        {
            trans_list_forward[j]     = SG_MALLOC(T_STATES, len);
            trans_list_forward_val[j] = SG_MALLOC(float64_t, len);
        }
        else
        {
            trans_list_forward[j]     = NULL;
            trans_list_forward_val[j] = NULL;
        }
    }

    for (int32_t i=0; i<num_trans; i++)
    {
        int32_t from = (int32_t)a_trans[i+num_trans] ;
        int32_t to   = (int32_t)a_trans[i] ;
        float64_t val = a_trans[i+num_trans*2] ;

        ASSERT(from>=0 && from<N)
        ASSERT(to>=0 && to<N)

        trans_list_forward[from][trans_list_forward_cnt[from]]=to ;
        trans_list_forward_val[from][trans_list_forward_cnt[from]]=val ;
        trans_list_forward_cnt[from]++ ;
        //ASSERT(trans_list_forward_cnt[from]<3000)
    } ;

    transition_matrix_a=NULL ;
    observation_matrix_b=NULL ;
    initial_state_distribution_p=p ;
    end_state_distribution_q=q ;
    transition_matrix_A=NULL ;
    observation_matrix_B=NULL ;

//  this->invalidate_model();
}


CHMM::CHMM(FILE* model_file, float64_t p_PSEUDO)
: CDistribution(), iterations(150), epsilon(1e-4), conv_it(5)
{
#ifdef USE_HMMPARALLEL_STRUCTURES
    SG_INFO("hmm is using %i separate tables\n",  parallel->get_num_threads())
#endif

    status=initialize_hmm(NULL, p_PSEUDO, model_file);
}

CHMM::~CHMM()
{
    SG_UNREF(p_observations);

    if (trans_list_forward_cnt)
      SG_FREE(trans_list_forward_cnt);
    if (trans_list_backward_cnt)
        SG_FREE(trans_list_backward_cnt);
    if (trans_list_forward)
    {
        for (int32_t i=0; i<trans_list_len; i++)
            if (trans_list_forward[i])
                SG_FREE(trans_list_forward[i]);
        SG_FREE(trans_list_forward);
    }
    if (trans_list_forward_val)
    {
        for (int32_t i=0; i<trans_list_len; i++)
            if (trans_list_forward_val[i])
                SG_FREE(trans_list_forward_val[i]);
        SG_FREE(trans_list_forward_val);
    }
    if (trans_list_backward)
      {
        for (int32_t i=0; i<trans_list_len; i++)
          if (trans_list_backward[i])
        SG_FREE(trans_list_backward[i]);
        SG_FREE(trans_list_backward);
      } ;

    free_state_dependend_arrays();

    if (!reused_caches)
    {
#ifdef USE_HMMPARALLEL_STRUCTURES
        if (mem_initialized)
        {
            for (int32_t i=0; i<parallel->get_num_threads(); i++)
            {
                SG_FREE(alpha_cache[i].table);
                SG_FREE(beta_cache[i].table);
                alpha_cache[i].table=NULL;
                beta_cache[i].table=NULL;
            }
        }
        SG_FREE(alpha_cache);
        SG_FREE(beta_cache);
        alpha_cache=NULL;
        beta_cache=NULL;
#else // USE_HMMPARALLEL_STRUCTURES
        SG_FREE(alpha_cache.table);
        SG_FREE(beta_cache.table);
        alpha_cache.table=NULL;
        beta_cache.table=NULL;
#endif // USE_HMMPARALLEL_STRUCTURES

        SG_FREE(states_per_observation_psi);
        states_per_observation_psi=NULL;
    }

#ifdef USE_LOGSUMARRAY
#ifdef USE_HMMPARALLEL_STRUCTURES
    {
        if (mem_initialized)
        {
            for (int32_t i=0; i<parallel->get_num_threads(); i++)
                SG_FREE(arrayS[i]);
        }
        SG_FREE(arrayS);
    } ;
#else //USE_HMMPARALLEL_STRUCTURES
    SG_FREE(arrayS);
#endif //USE_HMMPARALLEL_STRUCTURES
#endif //USE_LOGSUMARRAY

    if (!reused_caches)
    {
#ifdef USE_HMMPARALLEL_STRUCTURES
        if (mem_initialized)
        {
            SG_FREE(path_prob_updated);
            SG_FREE(path_prob_dimension);
            for (int32_t i=0; i<parallel->get_num_threads(); i++)
                SG_FREE(path[i]);
        }
#endif //USE_HMMPARALLEL_STRUCTURES
        SG_FREE(path);
    }
}

bool CHMM::train(CFeatures* data)
{
    if (data)
    {
        if (data->get_feature_class() != C_STRING ||
                data->get_feature_type() != F_WORD)
        {
            SG_ERROR("Expected features of class string type word\n")
        }
        set_observations((CStringFeatures<uint16_t>*) data);
    }
    return baum_welch_viterbi_train(BW_NORMAL);
}

bool CHMM::alloc_state_dependend_arrays()
{

    if (!transition_matrix_a && !observation_matrix_b &&
        !initial_state_distribution_p && !end_state_distribution_q)
    {
        transition_matrix_a=SG_MALLOC(float64_t, N*N);
        observation_matrix_b=SG_MALLOC(float64_t, N*M);
        initial_state_distribution_p=SG_MALLOC(float64_t, N);
        end_state_distribution_q=SG_MALLOC(float64_t, N);
        init_model_random();
        convert_to_log();
    }

#ifdef USE_HMMPARALLEL_STRUCTURES
    for (int32_t i=0; i<parallel->get_num_threads(); i++)
    {
        arrayN1[i]=SG_MALLOC(float64_t, N);
        arrayN2[i]=SG_MALLOC(float64_t, N);
    }
#else //USE_HMMPARALLEL_STRUCTURES
    arrayN1=SG_MALLOC(float64_t, N);
    arrayN2=SG_MALLOC(float64_t, N);
#endif //USE_HMMPARALLEL_STRUCTURES

#ifdef LOG_SUMARRAY
#ifdef USE_HMMPARALLEL_STRUCTURES
    for (int32_t i=0; i<parallel->get_num_threads(); i++)
        arrayS[i]=SG_MALLOC(float64_t, (int32_t)(this->N/2+1));
#else //USE_HMMPARALLEL_STRUCTURES
    arrayS=SG_MALLOC(float64_t, (int32_t)(this->N/2+1));
#endif //USE_HMMPARALLEL_STRUCTURES
#endif //LOG_SUMARRAY
    transition_matrix_A=SG_MALLOC(float64_t, this->N*this->N);
    observation_matrix_B=SG_MALLOC(float64_t, this->N*this->M);

    if (p_observations)
    {
#ifdef USE_HMMPARALLEL_STRUCTURES
        if (alpha_cache[0].table!=NULL)
#else //USE_HMMPARALLEL_STRUCTURES
        if (alpha_cache.table!=NULL)
#endif //USE_HMMPARALLEL_STRUCTURES
            set_observations(p_observations);
        else
            set_observation_nocache(p_observations);
        SG_UNREF(p_observations);
    }

    this->invalidate_model();

    return ((transition_matrix_A != NULL) && (observation_matrix_B != NULL) &&
            (transition_matrix_a != NULL) && (observation_matrix_b != NULL) &&
            (initial_state_distribution_p != NULL) &&
            (end_state_distribution_q != NULL));
}

void CHMM::free_state_dependend_arrays()
{
#ifdef USE_HMMPARALLEL_STRUCTURES
    if (arrayN1 && arrayN2)
    {
        for (int32_t i=0; i<parallel->get_num_threads(); i++)
        {
            SG_FREE(arrayN1[i]);
            SG_FREE(arrayN2[i]);

            arrayN1[i]=NULL;
            arrayN2[i]=NULL;
        }
    }
#endif
    SG_FREE(arrayN1);
    SG_FREE(arrayN2);
    arrayN1=NULL;
    arrayN2=NULL;

    if (observation_matrix_b)
    {
        SG_FREE(transition_matrix_A);
        SG_FREE(observation_matrix_B);
        SG_FREE(transition_matrix_a);
        SG_FREE(observation_matrix_b);
        SG_FREE(initial_state_distribution_p);
        SG_FREE(end_state_distribution_q);
    } ;

    transition_matrix_A=NULL;
    observation_matrix_B=NULL;
    transition_matrix_a=NULL;
    observation_matrix_b=NULL;
    initial_state_distribution_p=NULL;
    end_state_distribution_q=NULL;
}

bool CHMM::initialize_hmm(Model* m, float64_t pseudo, FILE* modelfile)
{
    //yes optimistic
    bool files_ok=true;

    trans_list_forward = NULL ;
    trans_list_forward_cnt = NULL ;
    trans_list_forward_val = NULL ;
    trans_list_backward = NULL ;
    trans_list_backward_cnt = NULL ;
    trans_list_len = 0 ;
    mem_initialized = false ;

    this->transition_matrix_a=NULL;
    this->observation_matrix_b=NULL;
    this->initial_state_distribution_p=NULL;
    this->end_state_distribution_q=NULL;
    this->PSEUDO= pseudo;
    this->model= m;
    this->p_observations=NULL;
    this->reused_caches=false;

#ifdef USE_HMMPARALLEL_STRUCTURES
    alpha_cache=SG_MALLOC(T_ALPHA_BETA, parallel->get_num_threads());
    beta_cache=SG_MALLOC(T_ALPHA_BETA, parallel->get_num_threads());
    states_per_observation_psi=SG_MALLOC(P_STATES, parallel->get_num_threads());

    for (int32_t i=0; i<parallel->get_num_threads(); i++)
    {
        this->alpha_cache[i].table=NULL;
        this->beta_cache[i].table=NULL;
        this->alpha_cache[i].dimension=0;
        this->beta_cache[i].dimension=0;
        this->states_per_observation_psi[i]=NULL ;
    }

#else // USE_HMMPARALLEL_STRUCTURES
    this->alpha_cache.table=NULL;
    this->beta_cache.table=NULL;
    this->alpha_cache.dimension=0;
    this->beta_cache.dimension=0;
    this->states_per_observation_psi=NULL ;
#endif //USE_HMMPARALLEL_STRUCTURES

    if (modelfile)
        files_ok= files_ok && load_model(modelfile);

#ifdef USE_HMMPARALLEL_STRUCTURES
    path_prob_updated=SG_MALLOC(bool, parallel->get_num_threads());
    path_prob_dimension=SG_MALLOC(int, parallel->get_num_threads());

    path=SG_MALLOC(P_STATES, parallel->get_num_threads());

    for (int32_t i=0; i<parallel->get_num_threads(); i++)
        this->path[i]=NULL;

#else // USE_HMMPARALLEL_STRUCTURES
    this->path=NULL;

#endif //USE_HMMPARALLEL_STRUCTURES

#ifdef USE_HMMPARALLEL_STRUCTURES
    arrayN1=SG_MALLOC(float64_t*, parallel->get_num_threads());
    arrayN2=SG_MALLOC(float64_t*, parallel->get_num_threads());
#endif //USE_HMMPARALLEL_STRUCTURES

#ifdef LOG_SUMARRAY
#ifdef USE_HMMPARALLEL_STRUCTURES
    arrayS=SG_MALLOC(float64_t*, parallel->get_num_threads());
#endif // USE_HMMPARALLEL_STRUCTURES
#endif //LOG_SUMARRAY

    alloc_state_dependend_arrays();

    this->loglikelihood=false;
    mem_initialized = true ;
    this->invalidate_model();

    return  ((files_ok) &&
            (transition_matrix_A != NULL) && (observation_matrix_B != NULL) &&
            (transition_matrix_a != NULL) && (observation_matrix_b != NULL) && (initial_state_distribution_p != NULL) &&
            (end_state_distribution_q != NULL));
}

//------------------------------------------------------------------------------------//

//forward algorithm
//calculates Pr[O_0,O_1, ..., O_t, q_time=S_i| lambda] for 0<= time <= T-1
//Pr[O|lambda] for time > T
float64_t CHMM::forward_comp(int32_t time, int32_t state, int32_t dimension)
{
    T_ALPHA_BETA_TABLE* alpha_new;
    T_ALPHA_BETA_TABLE* alpha;
    T_ALPHA_BETA_TABLE* dummy;
    if (time<1)
        time=0;


    int32_t wanted_time=time;

    if (ALPHA_CACHE(dimension).table)
    {
        alpha=&ALPHA_CACHE(dimension).table[0];
        alpha_new=&ALPHA_CACHE(dimension).table[N];
        time=p_observations->get_vector_length(dimension)+1;
    }
    else
    {
        alpha_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
        alpha=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
    }

    if (time<1)
        return get_p(state) + get_b(state, p_observations->get_feature(dimension,0));
    else
    {
        //initialization    alpha_1(i)=p_i*b_i(O_1)
        for (int32_t i=0; i<N; i++)
            alpha[i] = get_p(i) + get_b(i, p_observations->get_feature(dimension,0)) ;

        //induction     alpha_t+1(j) = (sum_i=1^N alpha_t(i)a_ij) b_j(O_t+1)
        for (register int32_t t=1; t<time && t < p_observations->get_vector_length(dimension); t++)
        {

            for (int32_t j=0; j<N; j++)
            {
                register int32_t i, num = trans_list_forward_cnt[j] ;
                float64_t sum=-CMath::INFTY;
                for (i=0; i < num; i++)
                {
                    int32_t ii = trans_list_forward[j][i] ;
                    sum = CMath::logarithmic_sum(sum, alpha[ii] + get_a(ii,j));
                } ;

                alpha_new[j]= sum + get_b(j, p_observations->get_feature(dimension,t));
            }

            if (!ALPHA_CACHE(dimension).table)
            {
                dummy=alpha;
                alpha=alpha_new;
                alpha_new=dummy;    //switch alpha/alpha_new
            }
            else
            {
                alpha=alpha_new;
                alpha_new+=N;       //perversely pointer arithmetic
            }
        }


        if (time<p_observations->get_vector_length(dimension))
        {
            register int32_t i, num=trans_list_forward_cnt[state];
            register float64_t sum=-CMath::INFTY;
            for (i=0; i<num; i++)
            {
                int32_t ii = trans_list_forward[state][i] ;
                sum= CMath::logarithmic_sum(sum, alpha[ii] + get_a(ii, state));
            } ;

            return sum + get_b(state, p_observations->get_feature(dimension,time));
        }
        else
        {
            // termination
            register int32_t i ;
            float64_t sum ;
            sum=-CMath::INFTY;
            for (i=0; i<N; i++)                     //sum over all paths
                sum=CMath::logarithmic_sum(sum, alpha[i] + get_q(i));   //to get model probability

            if (!ALPHA_CACHE(dimension).table)
                return sum;
            else
            {
                ALPHA_CACHE(dimension).dimension=dimension;
                ALPHA_CACHE(dimension).updated=true;
                ALPHA_CACHE(dimension).sum=sum;

                if (wanted_time<p_observations->get_vector_length(dimension))
                    return ALPHA_CACHE(dimension).table[wanted_time*N+state];
                else
                    return ALPHA_CACHE(dimension).sum;
            }
        }
    }
}


//forward algorithm
//calculates Pr[O_0,O_1, ..., O_t, q_time=S_i| lambda] for 0<= time <= T-1
//Pr[O|lambda] for time > T
float64_t CHMM::forward_comp_old(int32_t time, int32_t state, int32_t dimension)
{
    T_ALPHA_BETA_TABLE* alpha_new;
    T_ALPHA_BETA_TABLE* alpha;
    T_ALPHA_BETA_TABLE* dummy;
    if (time<1)
        time=0;

    int32_t wanted_time=time;

    if (ALPHA_CACHE(dimension).table)
    {
        alpha=&ALPHA_CACHE(dimension).table[0];
        alpha_new=&ALPHA_CACHE(dimension).table[N];
        time=p_observations->get_vector_length(dimension)+1;
    }
    else
    {
        alpha_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
        alpha=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
    }

    if (time<1)
        return get_p(state) + get_b(state, p_observations->get_feature(dimension,0));
    else
    {
        //initialization    alpha_1(i)=p_i*b_i(O_1)
        for (int32_t i=0; i<N; i++)
            alpha[i] = get_p(i) + get_b(i, p_observations->get_feature(dimension,0)) ;

        //induction     alpha_t+1(j) = (sum_i=1^N alpha_t(i)a_ij) b_j(O_t+1)
        for (register int32_t t=1; t<time && t < p_observations->get_vector_length(dimension); t++)
        {

            for (int32_t j=0; j<N; j++)
            {
                register int32_t i ;
#ifdef USE_LOGSUMARRAY
                for (i=0; i<(N>>1); i++)
                    ARRAYS(dimension)[i]=CMath::logarithmic_sum(alpha[i<<1] + get_a(i<<1,j),
                            alpha[(i<<1)+1] + get_a((i<<1)+1,j));
                if (N%2==1)
                    alpha_new[j]=get_b(j, p_observations->get_feature(dimension,t))+
                        CMath::logarithmic_sum(alpha[N-1]+get_a(N-1,j),
                                CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
                else
                    alpha_new[j]=get_b(j, p_observations->get_feature(dimension,t))+CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
#else //USE_LOGSUMARRAY
                float64_t sum=-CMath::INFTY;
                for (i=0; i<N; i++)
                    sum= CMath::logarithmic_sum(sum, alpha[i] + get_a(i,j));

                alpha_new[j]= sum + get_b(j, p_observations->get_feature(dimension,t));
#endif //USE_LOGSUMARRAY
            }

            if (!ALPHA_CACHE(dimension).table)
            {
                dummy=alpha;
                alpha=alpha_new;
                alpha_new=dummy;    //switch alpha/alpha_new
            }
            else
            {
                alpha=alpha_new;
                alpha_new+=N;       //perversely pointer arithmetic
            }
        }


        if (time<p_observations->get_vector_length(dimension))
        {
            register int32_t i;
#ifdef USE_LOGSUMARRAY
            for (i=0; i<(N>>1); i++)
                ARRAYS(dimension)[i]=CMath::logarithmic_sum(alpha[i<<1] + get_a(i<<1,state),
                        alpha[(i<<1)+1] + get_a((i<<1)+1,state));
            if (N%2==1)
                return get_b(state, p_observations->get_feature(dimension,time))+
                    CMath::logarithmic_sum(alpha[N-1]+get_a(N-1,state),
                            CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
            else
                return get_b(state, p_observations->get_feature(dimension,time))+CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
#else //USE_LOGSUMARRAY
            register float64_t sum=-CMath::INFTY;
            for (i=0; i<N; i++)
                sum= CMath::logarithmic_sum(sum, alpha[i] + get_a(i, state));

            return sum + get_b(state, p_observations->get_feature(dimension,time));
#endif //USE_LOGSUMARRAY
        }
        else
        {
            // termination
            register int32_t i ;
            float64_t sum ;
#ifdef USE_LOGSUMARRAY
            for (i=0; i<(N>>1); i++)
                ARRAYS(dimension)[i]=CMath::logarithmic_sum(alpha[i<<1] + get_q(i<<1),
                        alpha[(i<<1)+1] + get_q((i<<1)+1));
            if (N%2==1)
                sum=CMath::logarithmic_sum(alpha[N-1]+get_q(N-1),
                        CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
            else
                sum=CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
#else //USE_LOGSUMARRAY
            sum=-CMath::INFTY;
            for (i=0; i<N; i++)                               //sum over all paths
                sum=CMath::logarithmic_sum(sum, alpha[i] + get_q(i));     //to get model probability
#endif //USE_LOGSUMARRAY

            if (!ALPHA_CACHE(dimension).table)
                return sum;
            else
            {
                ALPHA_CACHE(dimension).dimension=dimension;
                ALPHA_CACHE(dimension).updated=true;
                ALPHA_CACHE(dimension).sum=sum;

                if (wanted_time<p_observations->get_vector_length(dimension))
                    return ALPHA_CACHE(dimension).table[wanted_time*N+state];
                else
                    return ALPHA_CACHE(dimension).sum;
            }
        }
    }
}


//backward algorithm
//calculates Pr[O_t+1,O_t+2, ..., O_T| q_time=S_i, lambda] for 0<= time <= T-1
//Pr[O|lambda] for time >= T
float64_t CHMM::backward_comp(int32_t time, int32_t state, int32_t dimension)
{
  T_ALPHA_BETA_TABLE* beta_new;
  T_ALPHA_BETA_TABLE* beta;
  T_ALPHA_BETA_TABLE* dummy;
  int32_t wanted_time=time;

  if (time<0)
    forward(time, state, dimension);

  if (BETA_CACHE(dimension).table)
    {
      beta=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-1)];
      beta_new=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-2)];
      time=-1;
    }
  else
    {
      beta_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
      beta=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
    }

  if (time>=p_observations->get_vector_length(dimension)-1)
    //    return 0;
    //  else if (time==p_observations->get_vector_length(dimension)-1)
    return get_q(state);
  else
    {
      //initialization  beta_T(i)=q(i)
      for (register int32_t i=0; i<N; i++)
    beta[i]=get_q(i);

      //induction       beta_t(i) = (sum_j=1^N a_ij*b_j(O_t+1)*beta_t+1(j)
      for (register int32_t t=p_observations->get_vector_length(dimension)-1; t>time+1 && t>0; t--)
    {
      for (register int32_t i=0; i<N; i++)
        {
          register int32_t j, num=trans_list_backward_cnt[i] ;
          float64_t sum=-CMath::INFTY;
          for (j=0; j<num; j++)
        {
          int32_t jj = trans_list_backward[i][j] ;
          sum= CMath::logarithmic_sum(sum, get_a(i, jj) + get_b(jj, p_observations->get_feature(dimension,t)) + beta[jj]);
        } ;
          beta_new[i]=sum;
        }

      if (!BETA_CACHE(dimension).table)
        {
          dummy=beta;
          beta=beta_new;
          beta_new=dummy;   //switch beta/beta_new
        }
      else
        {
          beta=beta_new;
          beta_new-=N;      //perversely pointer arithmetic
        }
    }

      if (time>=0)
    {
      register int32_t j, num=trans_list_backward_cnt[state] ;
      float64_t sum=-CMath::INFTY;
      for (j=0; j<num; j++)
        {
          int32_t jj = trans_list_backward[state][j] ;
          sum= CMath::logarithmic_sum(sum, get_a(state, jj) + get_b(jj, p_observations->get_feature(dimension,time+1))+beta[jj]);
        } ;
      return sum;
    }
      else // time<0
    {
      if (BETA_CACHE(dimension).table)
        {
          float64_t sum=-CMath::INFTY;
          for (register int32_t j=0; j<N; j++)
        sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
          BETA_CACHE(dimension).sum=sum;
          BETA_CACHE(dimension).dimension=dimension;
          BETA_CACHE(dimension).updated=true;

          if (wanted_time<p_observations->get_vector_length(dimension))
        return BETA_CACHE(dimension).table[wanted_time*N+state];
          else
        return BETA_CACHE(dimension).sum;
        }
      else
        {
          float64_t sum=-CMath::INFTY; // apply LOG_SUM_ARRAY -- no cache ... does not make very sense anyway...
          for (register int32_t j=0; j<N; j++)
        sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
          return sum;
        }
    }
    }
}


float64_t CHMM::backward_comp_old(
    int32_t time, int32_t state, int32_t dimension)
{
    T_ALPHA_BETA_TABLE* beta_new;
    T_ALPHA_BETA_TABLE* beta;
    T_ALPHA_BETA_TABLE* dummy;
    int32_t wanted_time=time;

    if (time<0)
        forward(time, state, dimension);

    if (BETA_CACHE(dimension).table)
    {
        beta=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-1)];
        beta_new=&BETA_CACHE(dimension).table[N*(p_observations->get_vector_length(dimension)-2)];
        time=-1;
    }
    else
    {
        beta_new=(T_ALPHA_BETA_TABLE*)ARRAYN1(dimension);
        beta=(T_ALPHA_BETA_TABLE*)ARRAYN2(dimension);
    }

    if (time>=p_observations->get_vector_length(dimension)-1)
        //    return 0;
        //  else if (time==p_observations->get_vector_length(dimension)-1)
        return get_q(state);
    else
    {
        //initialization    beta_T(i)=q(i)
        for (register int32_t i=0; i<N; i++)
            beta[i]=get_q(i);

        //induction     beta_t(i) = (sum_j=1^N a_ij*b_j(O_t+1)*beta_t+1(j)
        for (register int32_t t=p_observations->get_vector_length(dimension)-1; t>time+1 && t>0; t--)
        {
            for (register int32_t i=0; i<N; i++)
            {
                register int32_t j ;
#ifdef USE_LOGSUMARRAY
                for (j=0; j<(N>>1); j++)
                    ARRAYS(dimension)[j]=CMath::logarithmic_sum(
                            get_a(i, j<<1) + get_b(j<<1, p_observations->get_feature(dimension,t)) + beta[j<<1],
                            get_a(i, (j<<1)+1) + get_b((j<<1)+1, p_observations->get_feature(dimension,t)) + beta[(j<<1)+1]);
                if (N%2==1)
                    beta_new[i]=CMath::logarithmic_sum(get_a(i, N-1) + get_b(N-1, p_observations->get_feature(dimension,t)) + beta[N-1],
                            CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
                else
                    beta_new[i]=CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
#else //USE_LOGSUMARRAY
                float64_t sum=-CMath::INFTY;
                for (j=0; j<N; j++)
                    sum= CMath::logarithmic_sum(sum, get_a(i, j) + get_b(j, p_observations->get_feature(dimension,t)) + beta[j]);

                beta_new[i]=sum;
#endif //USE_LOGSUMARRAY
            }

            if (!BETA_CACHE(dimension).table)
            {
                dummy=beta;
                beta=beta_new;
                beta_new=dummy; //switch beta/beta_new
            }
            else
            {
                beta=beta_new;
                beta_new-=N;        //perversely pointer arithmetic
            }
        }

        if (time>=0)
        {
            register int32_t j ;
#ifdef USE_LOGSUMARRAY
            for (j=0; j<(N>>1); j++)
                ARRAYS(dimension)[j]=CMath::logarithmic_sum(
                        get_a(state, j<<1) + get_b(j<<1, p_observations->get_feature(dimension,time+1)) + beta[j<<1],
                        get_a(state, (j<<1)+1) + get_b((j<<1)+1, p_observations->get_feature(dimension,time+1)) + beta[(j<<1)+1]);
            if (N%2==1)
                return CMath::logarithmic_sum(get_a(state, N-1) + get_b(N-1, p_observations->get_feature(dimension,time+1)) + beta[N-1],
                        CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
            else
                return CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
#else //USE_LOGSUMARRAY
            float64_t sum=-CMath::INFTY;
            for (j=0; j<N; j++)
                sum= CMath::logarithmic_sum(sum, get_a(state, j) + get_b(j, p_observations->get_feature(dimension,time+1))+beta[j]);

            return sum;
#endif //USE_LOGSUMARRAY
        }
        else // time<0
        {
            if (BETA_CACHE(dimension).table)
            {
#ifdef USE_LOGSUMARRAY//AAA
                for (int32_t j=0; j<(N>>1); j++)
                    ARRAYS(dimension)[j]=CMath::logarithmic_sum(get_p(j<<1) + get_b(j<<1, p_observations->get_feature(dimension,0))+beta[j<<1],
                            get_p((j<<1)+1) + get_b((j<<1)+1, p_observations->get_feature(dimension,0))+beta[(j<<1)+1]) ;
                if (N%2==1)
                    BETA_CACHE(dimension).sum=CMath::logarithmic_sum(get_p(N-1) + get_b(N-1, p_observations->get_feature(dimension,0))+beta[N-1],
                            CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1)) ;
                else
                    BETA_CACHE(dimension).sum=CMath::logarithmic_sum_array(ARRAYS(dimension), N>>1) ;
#else //USE_LOGSUMARRAY
                float64_t sum=-CMath::INFTY;
                for (register int32_t j=0; j<N; j++)
                    sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
                BETA_CACHE(dimension).sum=sum;
#endif //USE_LOGSUMARRAY
                BETA_CACHE(dimension).dimension=dimension;
                BETA_CACHE(dimension).updated=true;

                if (wanted_time<p_observations->get_vector_length(dimension))
                    return BETA_CACHE(dimension).table[wanted_time*N+state];
                else
                    return BETA_CACHE(dimension).sum;
            }
            else
            {
                float64_t sum=-CMath::INFTY; // apply LOG_SUM_ARRAY -- no cache ... does not make very sense anyway...
                for (register int32_t j=0; j<N; j++)
                    sum= CMath::logarithmic_sum(sum, get_p(j) + get_b(j, p_observations->get_feature(dimension,0))+beta[j]);
                return sum;
            }
        }
    }
}

//calculates probability  of best path through the model lambda AND path itself
//using viterbi algorithm
float64_t CHMM::best_path(int32_t dimension)
{
    if (!p_observations)
        return -1;

    if (dimension==-1)
    {
        if (!all_path_prob_updated)
        {
            SG_INFO("computing full viterbi likelihood\n")
            float64_t sum = 0 ;
            for (int32_t i=0; i<p_observations->get_num_vectors(); i++)
                sum+=best_path(i) ;
            sum /= p_observations->get_num_vectors() ;
            all_pat_prob=sum ;
            all_path_prob_updated=true ;
            return sum ;
        } else
            return all_pat_prob ;
    } ;

    if (!STATES_PER_OBSERVATION_PSI(dimension))
        return -1 ;

    if (dimension >= p_observations->get_num_vectors())
        return -1;

    if (PATH_PROB_UPDATED(dimension) && dimension==PATH_PROB_DIMENSION(dimension))
        return pat_prob;
    else
    {
        register float64_t* delta= ARRAYN2(dimension);
        register float64_t* delta_new= ARRAYN1(dimension);

        { //initialization
            for (register int32_t i=0; i<N; i++)
            {
                delta[i]=get_p(i)+get_b(i, p_observations->get_feature(dimension,0));
                set_psi(0, i, 0, dimension);
            }
        }

#ifdef USE_PATHDEBUG
        float64_t worst=-CMath::INFTY/4 ;
#endif
        //recursion
        for (register int32_t t=1; t<p_observations->get_vector_length(dimension); t++)
        {
            register float64_t* dummy;
            register int32_t NN=N ;
            for (register int32_t j=0; j<NN; j++)
            {
                register float64_t * matrix_a=&transition_matrix_a[j*N] ; // sorry for that
                register float64_t maxj=delta[0] + matrix_a[0];
                register int32_t argmax=0;

                for (register int32_t i=1; i<NN; i++)
                {
                    register float64_t temp = delta[i] + matrix_a[i];

                    if (temp>maxj)
                    {
                        maxj=temp;
                        argmax=i;
                    }
                }
#ifdef FIX_POS
                if ((!model) || (model->get_fix_pos_state(t,j,NN)!=Model::FIX_DISALLOWED))
#endif
                    delta_new[j]=maxj + get_b(j,p_observations->get_feature(dimension,t));
#ifdef FIX_POS
                else
                    delta_new[j]=maxj + get_b(j,p_observations->get_feature(dimension,t)) + Model::DISALLOWED_PENALTY;
#endif
                set_psi(t, j, argmax, dimension);
            }

#ifdef USE_PATHDEBUG
            float64_t best=log(0) ;
            for (int32_t jj=0; jj<N; jj++)
                if (delta_new[jj]>best)
                    best=delta_new[jj] ;

            if (best<-CMath::INFTY/2)
            {
                SG_DEBUG("worst case at %i: %e:%e\n", t, best, worst)
                worst=best ;
            } ;
#endif

            dummy=delta;
            delta=delta_new;
            delta_new=dummy;    //switch delta/delta_new
        }


        { //termination
            register float64_t maxj=delta[0]+get_q(0);
            register int32_t argmax=0;

            for (register int32_t i=1; i<N; i++)
            {
                register float64_t temp=delta[i]+get_q(i);

                if (temp>maxj)
                {
                    maxj=temp;
                    argmax=i;
                }
            }
            pat_prob=maxj;
            PATH(dimension)[p_observations->get_vector_length(dimension)-1]=argmax;
        } ;


        { //state sequence backtracking
            for (register int32_t t=p_observations->get_vector_length(dimension)-1; t>0; t--)
            {
                PATH(dimension)[t-1]=get_psi(t, PATH(dimension)[t], dimension);
            }
        }
        PATH_PROB_UPDATED(dimension)=true;
        PATH_PROB_DIMENSION(dimension)=dimension;
        return pat_prob ;
    }
}

#ifndef USE_HMMPARALLEL
float64_t CHMM::model_probability_comp()
{
    //for faster calculation cache model probability
    mod_prob=0 ;
    for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++) //sum in log space
        mod_prob+=forward(p_observations->get_vector_length(dim), 0, dim);

    mod_prob_updated=true;
    return mod_prob;
}

#else

float64_t CHMM::model_probability_comp()
{
    pthread_t *threads=SG_MALLOC(pthread_t, parallel->get_num_threads());
    S_BW_THREAD_PARAM *params=SG_MALLOC(S_BW_THREAD_PARAM, parallel->get_num_threads());

    SG_INFO("computing full model probablity\n")
    mod_prob=0;

    for (int32_t cpu=0; cpu<parallel->get_num_threads(); cpu++)
    {
        params[cpu].hmm=this ;
        params[cpu].dim_start= p_observations->get_num_vectors()*cpu/parallel->get_num_threads();
        params[cpu].dim_stop= p_observations->get_num_vectors()*(cpu+1)/parallel->get_num_threads();
        params[cpu].p_buf=SG_MALLOC(float64_t, N);
        params[cpu].q_buf=SG_MALLOC(float64_t, N);
        params[cpu].a_buf=SG_MALLOC(float64_t, N*N);
        params[cpu].b_buf=SG_MALLOC(float64_t, N*M);
        pthread_create(&threads[cpu], NULL, bw_dim_prefetch, (void*)&params[cpu]);
    }

    for (int32_t cpu=0; cpu<parallel->get_num_threads(); cpu++)
    {
        pthread_join(threads[cpu], NULL);
        mod_prob+=params[cpu].ret;
    }

    for (int32_t i=0; i<parallel->get_num_threads(); i++)
    {
        SG_FREE(params[i].p_buf);
        SG_FREE(params[i].q_buf);
        SG_FREE(params[i].a_buf);
        SG_FREE(params[i].b_buf);
    }

    SG_FREE(threads);
    SG_FREE(params);

    mod_prob_updated=true;
    return mod_prob;
}

void* CHMM::bw_dim_prefetch(void* params)
{
    CHMM* hmm=((S_BW_THREAD_PARAM*) params)->hmm;
    int32_t start=((S_BW_THREAD_PARAM*) params)->dim_start;
    int32_t stop=((S_BW_THREAD_PARAM*) params)->dim_stop;
    float64_t* p_buf=((S_BW_THREAD_PARAM*) params)->p_buf;
    float64_t* q_buf=((S_BW_THREAD_PARAM*) params)->q_buf;
    float64_t* a_buf=((S_BW_THREAD_PARAM*) params)->a_buf;
    float64_t* b_buf=((S_BW_THREAD_PARAM*) params)->b_buf;
    ((S_BW_THREAD_PARAM*)params)->ret=0;

    for (int32_t dim=start; dim<stop; dim++)
    {
        hmm->forward_comp(hmm->p_observations->get_vector_length(dim), hmm->N-1, dim) ;
        hmm->backward_comp(hmm->p_observations->get_vector_length(dim), hmm->N-1, dim) ;
        float64_t modprob=hmm->model_probability(dim) ;
        hmm->ab_buf_comp(p_buf, q_buf, a_buf, b_buf, dim) ;
        ((S_BW_THREAD_PARAM*)params)->ret+= modprob;
    }
    return NULL ;
}

void* CHMM::bw_single_dim_prefetch(void * params)
{
    CHMM* hmm=((S_BW_THREAD_PARAM*)params)->hmm ;
    int32_t dim=((S_DIM_THREAD_PARAM*)params)->dim ;
    ((S_DIM_THREAD_PARAM*)params)->prob_sum = hmm->model_probability(dim);
    return NULL ;
}

void* CHMM::vit_dim_prefetch(void * params)
{
    CHMM* hmm=((S_DIM_THREAD_PARAM*)params)->hmm ;
    int32_t dim=((S_DIM_THREAD_PARAM*)params)->dim ;
    ((S_DIM_THREAD_PARAM*)params)->prob_sum = hmm->best_path(dim);
    return NULL ;
}

#endif //USE_HMMPARALLEL


#ifdef USE_HMMPARALLEL

void CHMM::ab_buf_comp(
    float64_t* p_buf, float64_t* q_buf, float64_t *a_buf, float64_t* b_buf,
    int32_t dim)
{
    int32_t i,j,t ;
    float64_t a_sum;
    float64_t b_sum;

    float64_t dimmodprob=model_probability(dim);

    for (i=0; i<N; i++)
    {
        //estimate initial+end state distribution numerator
        p_buf[i]=get_p(i)+get_b(i,p_observations->get_feature(dim,0))+backward(0,i,dim) - dimmodprob;
        q_buf[i]=forward(p_observations->get_vector_length(dim)-1, i, dim)+get_q(i) - dimmodprob;

        //estimate a
        for (j=0; j<N; j++)
        {
            a_sum=-CMath::INFTY;

            for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
            {
                a_sum= CMath::logarithmic_sum(a_sum, forward(t,i,dim)+
                        get_a(i,j)+get_b(j,p_observations->get_feature(dim,t+1))+backward(t+1,j,dim));
            }
            a_buf[N*i+j]=a_sum-dimmodprob ;
        }

        //estimate b
        for (j=0; j<M; j++)
        {
            b_sum=-CMath::INFTY;

            for (t=0; t<p_observations->get_vector_length(dim); t++)
            {
                if (p_observations->get_feature(dim,t)==j)
                    b_sum=CMath::logarithmic_sum(b_sum, forward(t,i,dim)+backward(t, i, dim));
            }

            b_buf[M*i+j]=b_sum-dimmodprob ;
        }
    }
}

//estimates new model lambda out of lambda_train using baum welch algorithm
void CHMM::estimate_model_baum_welch(CHMM* hmm)
{
    int32_t i,j,cpu;
    float64_t fullmodprob=0;    //for all dims

    //clear actual model a,b,p,q are used as numerator
    for (i=0; i<N; i++)
    {
      if (hmm->get_p(i)>CMath::ALMOST_NEG_INFTY)
        set_p(i,log(PSEUDO));
      else
        set_p(i,hmm->get_p(i));
      if (hmm->get_q(i)>CMath::ALMOST_NEG_INFTY)
        set_q(i,log(PSEUDO));
      else
        set_q(i,hmm->get_q(i));

      for (j=0; j<N; j++)
        if (hmm->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
          set_a(i,j, log(PSEUDO));
        else
          set_a(i,j,hmm->get_a(i,j));
      for (j=0; j<M; j++)
        if (hmm->get_b(i,j)>CMath::ALMOST_NEG_INFTY)
          set_b(i,j, log(PSEUDO));
        else
          set_b(i,j,hmm->get_b(i,j));
    }
    invalidate_model();

    int32_t num_threads = parallel->get_num_threads();

    pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
    S_BW_THREAD_PARAM *params=SG_MALLOC(S_BW_THREAD_PARAM, num_threads);

    if (p_observations->get_num_vectors()<num_threads)
        num_threads=p_observations->get_num_vectors();

    for (cpu=0; cpu<num_threads; cpu++)
    {
        params[cpu].p_buf=SG_MALLOC(float64_t, N);
        params[cpu].q_buf=SG_MALLOC(float64_t, N);
        params[cpu].a_buf=SG_MALLOC(float64_t, N*N);
        params[cpu].b_buf=SG_MALLOC(float64_t, N*M);

        params[cpu].hmm=hmm;
        int32_t start = p_observations->get_num_vectors()*cpu / num_threads;
        int32_t stop=p_observations->get_num_vectors()*(cpu+1) / num_threads;

        if (cpu == parallel->get_num_threads()-1)
            stop=p_observations->get_num_vectors();

        ASSERT(start<stop)
        params[cpu].dim_start=start;
        params[cpu].dim_stop=stop;

        pthread_create(&threads[cpu], NULL, bw_dim_prefetch, &params[cpu]);
    }

    for (cpu=0; cpu<num_threads; cpu++)
    {
        pthread_join(threads[cpu], NULL);

        for (i=0; i<N; i++)
        {
            //estimate initial+end state distribution numerator
            set_p(i, CMath::logarithmic_sum(get_p(i), params[cpu].p_buf[i]));
            set_q(i, CMath::logarithmic_sum(get_q(i), params[cpu].q_buf[i]));

            //estimate numerator for a
            for (j=0; j<N; j++)
                set_a(i,j, CMath::logarithmic_sum(get_a(i,j), params[cpu].a_buf[N*i+j]));

            //estimate numerator for b
            for (j=0; j<M; j++)
                set_b(i,j, CMath::logarithmic_sum(get_b(i,j), params[cpu].b_buf[M*i+j]));
        }

        fullmodprob+=params[cpu].ret;

    }

    for (cpu=0; cpu<num_threads; cpu++)
    {
        SG_FREE(params[cpu].p_buf);
        SG_FREE(params[cpu].q_buf);
        SG_FREE(params[cpu].a_buf);
        SG_FREE(params[cpu].b_buf);
    }

    SG_FREE(threads);
    SG_FREE(params);

    //cache hmm model probability
    hmm->mod_prob=fullmodprob;
    hmm->mod_prob_updated=true ;

    //new model probability is unknown
    normalize();
    invalidate_model();
}

#else // USE_HMMPARALLEL

//estimates new model lambda out of lambda_estimate using baum welch algorithm
void CHMM::estimate_model_baum_welch(CHMM* estimate)
{
    int32_t i,j,t,dim;
    float64_t a_sum, b_sum; //numerator
    float64_t dimmodprob=0; //model probability for dim
    float64_t fullmodprob=0;    //for all dims

    //clear actual model a,b,p,q are used as numerator
    for (i=0; i<N; i++)
    {
        if (estimate->get_p(i)>CMath::ALMOST_NEG_INFTY)
            set_p(i,log(PSEUDO));
        else
            set_p(i,estimate->get_p(i));
        if (estimate->get_q(i)>CMath::ALMOST_NEG_INFTY)
            set_q(i,log(PSEUDO));
        else
            set_q(i,estimate->get_q(i));

        for (j=0; j<N; j++)
            if (estimate->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
                set_a(i,j, log(PSEUDO));
            else
                set_a(i,j,estimate->get_a(i,j));
        for (j=0; j<M; j++)
            if (estimate->get_b(i,j)>CMath::ALMOST_NEG_INFTY)
                set_b(i,j, log(PSEUDO));
            else
                set_b(i,j,estimate->get_b(i,j));
    }
    invalidate_model();

    //change summation order to make use of alpha/beta caches
    for (dim=0; dim<p_observations->get_num_vectors(); dim++)
    {
        dimmodprob=estimate->model_probability(dim);
        fullmodprob+=dimmodprob ;

        for (i=0; i<N; i++)
        {
            //estimate initial+end state distribution numerator
            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));
            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 ));

            int32_t num = trans_list_backward_cnt[i] ;

            //estimate a
            for (j=0; j<num; j++)
            {
                int32_t jj = trans_list_backward[i][j] ;
                a_sum=-CMath::INFTY;

                for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
                {
                    a_sum= CMath::logarithmic_sum(a_sum, estimate->forward(t,i,dim)+
                            estimate->get_a(i,jj)+estimate->get_b(jj,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,jj,dim));
                }
                set_a(i,jj, CMath::logarithmic_sum(get_a(i,jj), a_sum-dimmodprob));
            }

            //estimate b
            for (j=0; j<M; j++)
            {
                b_sum=-CMath::INFTY;

                for (t=0; t<p_observations->get_vector_length(dim); t++)
                {
                    if (p_observations->get_feature(dim,t)==j)
                        b_sum=CMath::logarithmic_sum(b_sum, estimate->forward(t,i,dim)+estimate->backward(t, i, dim));
                }

                set_b(i,j, CMath::logarithmic_sum(get_b(i,j), b_sum-dimmodprob));
            }
        }
    }

    //cache estimate model probability
    estimate->mod_prob=fullmodprob;
    estimate->mod_prob_updated=true ;

    //new model probability is unknown
    normalize();
    invalidate_model();
}

//estimates new model lambda out of lambda_estimate using baum welch algorithm
void CHMM::estimate_model_baum_welch_old(CHMM* estimate)
{
    int32_t i,j,t,dim;
    float64_t a_sum, b_sum; //numerator
    float64_t dimmodprob=0; //model probability for dim
    float64_t fullmodprob=0;    //for all dims

    //clear actual model a,b,p,q are used as numerator
    for (i=0; i<N; i++)
      {
        if (estimate->get_p(i)>CMath::ALMOST_NEG_INFTY)
          set_p(i,log(PSEUDO));
        else
          set_p(i,estimate->get_p(i));
        if (estimate->get_q(i)>CMath::ALMOST_NEG_INFTY)
          set_q(i,log(PSEUDO));
        else
          set_q(i,estimate->get_q(i));

        for (j=0; j<N; j++)
          if (estimate->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
        set_a(i,j, log(PSEUDO));
          else
        set_a(i,j,estimate->get_a(i,j));
        for (j=0; j<M; j++)
          if (estimate->get_b(i,j)>CMath::ALMOST_NEG_INFTY)
        set_b(i,j, log(PSEUDO));
          else
        set_b(i,j,estimate->get_b(i,j));
      }
    invalidate_model();

    //change summation order to make use of alpha/beta caches
    for (dim=0; dim<p_observations->get_num_vectors(); dim++)
      {
        dimmodprob=estimate->model_probability(dim);
        fullmodprob+=dimmodprob ;

        for (i=0; i<N; i++)
          {
        //estimate initial+end state distribution numerator
        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));
        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 ));

        //estimate a
        for (j=0; j<N; j++)
          {
            a_sum=-CMath::INFTY;

            for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
              {
            a_sum= CMath::logarithmic_sum(a_sum, estimate->forward(t,i,dim)+
                            estimate->get_a(i,j)+estimate->get_b(j,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,j,dim));
              }
            set_a(i,j, CMath::logarithmic_sum(get_a(i,j), a_sum-dimmodprob));
          }

        //estimate b
        for (j=0; j<M; j++)
          {
            b_sum=-CMath::INFTY;

            for (t=0; t<p_observations->get_vector_length(dim); t++)
              {
            if (p_observations->get_feature(dim,t)==j)
              b_sum=CMath::logarithmic_sum(b_sum, estimate->forward(t,i,dim)+estimate->backward(t, i, dim));
              }

            set_b(i,j, CMath::logarithmic_sum(get_b(i,j), b_sum-dimmodprob));
          }
          }
      }

    //cache estimate model probability
    estimate->mod_prob=fullmodprob;
    estimate->mod_prob_updated=true ;

    //new model probability is unknown
    normalize();
    invalidate_model();
}
#endif // USE_HMMPARALLEL

//estimates new model lambda out of lambda_estimate using baum welch algorithm
// optimize only p, q, a but not b
void CHMM::estimate_model_baum_welch_trans(CHMM* estimate)
{
    int32_t i,j,t,dim;
    float64_t a_sum;    //numerator
    float64_t dimmodprob=0; //model probability for dim
    float64_t fullmodprob=0;    //for all dims

    //clear actual model a,b,p,q are used as numerator
    for (i=0; i<N; i++)
      {
        if (estimate->get_p(i)>CMath::ALMOST_NEG_INFTY)
          set_p(i,log(PSEUDO));
        else
          set_p(i,estimate->get_p(i));
        if (estimate->get_q(i)>CMath::ALMOST_NEG_INFTY)
          set_q(i,log(PSEUDO));
        else
          set_q(i,estimate->get_q(i));

        for (j=0; j<N; j++)
          if (estimate->get_a(i,j)>CMath::ALMOST_NEG_INFTY)
        set_a(i,j, log(PSEUDO));
          else
        set_a(i,j,estimate->get_a(i,j));
        for (j=0; j<M; j++)
          set_b(i,j,estimate->get_b(i,j));
      }
    invalidate_model();

    //change summation order to make use of alpha/beta caches
    for (dim=0; dim<p_observations->get_num_vectors(); dim++)
      {
        dimmodprob=estimate->model_probability(dim);
        fullmodprob+=dimmodprob ;

        for (i=0; i<N; i++)
          {
        //estimate initial+end state distribution numerator
        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));
        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 ));

        int32_t num = trans_list_backward_cnt[i] ;
        //estimate a
        for (j=0; j<num; j++)
          {
            int32_t jj = trans_list_backward[i][j] ;
            a_sum=-CMath::INFTY;

            for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
              {
            a_sum= CMath::logarithmic_sum(a_sum, estimate->forward(t,i,dim)+
                            estimate->get_a(i,jj)+estimate->get_b(jj,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,jj,dim));
              }
            set_a(i,jj, CMath::logarithmic_sum(get_a(i,jj), a_sum-dimmodprob));
          }
          }
      }

    //cache estimate model probability
    estimate->mod_prob=fullmodprob;
    estimate->mod_prob_updated=true ;

    //new model probability is unknown
    normalize();
    invalidate_model();
}



//estimates new model lambda out of lambda_estimate using baum welch algorithm
void CHMM::estimate_model_baum_welch_defined(CHMM* estimate)
{
    int32_t i,j,old_i,k,t,dim;
    float64_t a_sum_num, b_sum_num;     //numerator
    float64_t a_sum_denom, b_sum_denom; //denominator
    float64_t dimmodprob=-CMath::INFTY; //model probability for dim
    float64_t fullmodprob=0;            //for all dims
    float64_t* A=ARRAYN1(0);
    float64_t* B=ARRAYN2(0);

    //clear actual model a,b,p,q are used as numerator
    //A,B as denominator for a,b
    for (k=0; (i=model->get_learn_p(k))!=-1; k++)
        set_p(i,log(PSEUDO));

    for (k=0; (i=model->get_learn_q(k))!=-1; k++)
        set_q(i,log(PSEUDO));

    for (k=0; (i=model->get_learn_a(k,0))!=-1; k++)
    {
        j=model->get_learn_a(k,1);
        set_a(i,j, log(PSEUDO));
    }

    for (k=0; (i=model->get_learn_b(k,0))!=-1; k++)
    {
        j=model->get_learn_b(k,1);
        set_b(i,j, log(PSEUDO));
    }

    for (i=0; i<N; i++)
    {
        A[i]=log(PSEUDO);
        B[i]=log(PSEUDO);
    }

#ifdef USE_HMMPARALLEL
    int32_t num_threads = parallel->get_num_threads();
    pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
    S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);

    if (p_observations->get_num_vectors()<num_threads)
        num_threads=p_observations->get_num_vectors();
#endif

    //change summation order to make use of alpha/beta caches
    for (dim=0; dim<p_observations->get_num_vectors(); dim++)
    {
#ifdef USE_HMMPARALLEL
        if (dim%num_threads==0)
        {
            for (i=0; i<num_threads; i++)
            {
                if (dim+i<p_observations->get_num_vectors())
                {
                    params[i].hmm=estimate ;
                    params[i].dim=dim+i ;
                    pthread_create(&threads[i], NULL, bw_single_dim_prefetch, (void*)&params[i]) ;
                }
            }
            for (i=0; i<num_threads; i++)
            {
                if (dim+i<p_observations->get_num_vectors())
                {
                    pthread_join(threads[i], NULL);
                    dimmodprob = params[i].prob_sum;
                }
            }
        }
#else
        dimmodprob=estimate->model_probability(dim);
#endif // USE_HMMPARALLEL

        //and denominator
        fullmodprob+= dimmodprob;

        //estimate initial+end state distribution numerator
        for (k=0; (i=model->get_learn_p(k))!=-1; k++)
            set_p(i, CMath::logarithmic_sum(get_p(i), estimate->forward(0,i,dim)+estimate->backward(0,i,dim) - dimmodprob ) );

        for (k=0; (i=model->get_learn_q(k))!=-1; k++)
            set_q(i, CMath::logarithmic_sum(get_q(i), estimate->forward(p_observations->get_vector_length(dim)-1, i, dim)+
                        estimate->backward(p_observations->get_vector_length(dim)-1, i, dim)  - dimmodprob ) );

        //estimate a
        old_i=-1;
        for (k=0; (i=model->get_learn_a(k,0))!=-1; k++)
        {
            //denominator is constant for j
            //therefore calculate it first
            if (old_i!=i)
            {
                old_i=i;
                a_sum_denom=-CMath::INFTY;

                for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
                    a_sum_denom= CMath::logarithmic_sum(a_sum_denom, estimate->forward(t,i,dim)+estimate->backward(t,i,dim));

                A[i]= CMath::logarithmic_sum(A[i], a_sum_denom-dimmodprob);
            }

            j=model->get_learn_a(k,1);
            a_sum_num=-CMath::INFTY;
            for (t=0; t<p_observations->get_vector_length(dim)-1; t++)
            {
                a_sum_num= CMath::logarithmic_sum(a_sum_num, estimate->forward(t,i,dim)+
                        estimate->get_a(i,j)+estimate->get_b(j,p_observations->get_feature(dim,t+1))+estimate->backward(t+1,j,dim));
            }

            set_a(i,j, CMath::logarithmic_sum(get_a(i,j), a_sum_num-dimmodprob));
        }

        //estimate  b
        old_i=-1;
        for (k=0; (i=model->get_learn_b(k,0))!=-1; k++)
        {

            //denominator is constant for j
            //therefore calculate it first
            if (old_i!=i)
            {
                old_i=i;
                b_sum_denom=-CMath::INFTY;

                for (t=0; t<p_observations->get_vector_length(dim); t++)
                    b_sum_denom= CMath::logarithmic_sum(b_sum_denom, estimate->forward(t,i,dim)+estimate->backward(t,i,dim));

                B[i]= CMath::logarithmic_sum(B[i], b_sum_denom-dimmodprob);
            }

            j=model->get_learn_b(k,1);
            b_sum_num=-CMath::INFTY;
            for (t=0; t<p_observations->get_vector_length(dim); t++)
            {
                if (p_observations->get_feature(dim,t)==j)
                    b_sum_num=CMath::logarithmic_sum(b_sum_num, estimate->forward(t,i,dim)+estimate->backward(t, i, dim));
            }

            set_b(i,j, CMath::logarithmic_sum(get_b(i,j), b_sum_num-dimmodprob));
        }
    }
#ifdef USE_HMMPARALLEL
    SG_FREE(threads);
    SG_FREE(params);
#endif


    //calculate estimates
    for (k=0; (i=model->get_learn_p(k))!=-1; k++)
        set_p(i, get_p(i)-log(p_observations->get_num_vectors()+N*PSEUDO) );

    for (k=0; (i=model->get_learn_q(k))!=-1; k++)
        set_q(i, get_q(i)-log(p_observations->get_num_vectors()+N*PSEUDO) );

    for (k=0; (i=model->get_learn_a(k,0))!=-1; k++)
    {
        j=model->get_learn_a(k,1);
        set_a(i,j, get_a(i,j) - A[i]);
    }

    for (k=0; (i=model->get_learn_b(k,0))!=-1; k++)
    {
        j=model->get_learn_b(k,1);
        set_b(i,j, get_b(i,j) - B[i]);
    }

    //cache estimate model probability
    estimate->mod_prob=fullmodprob;
    estimate->mod_prob_updated=true ;

    //new model probability is unknown
    normalize();
    invalidate_model();
}

//estimates new model lambda out of lambda_estimate using viterbi algorithm
void CHMM::estimate_model_viterbi(CHMM* estimate)
{
    int32_t i,j,t;
    float64_t sum;
    float64_t* P=ARRAYN1(0);
    float64_t* Q=ARRAYN2(0);

    path_deriv_updated=false ;

    //initialize with pseudocounts
    for (i=0; i<N; i++)
    {
        for (j=0; j<N; j++)
            set_A(i,j, PSEUDO);

        for (j=0; j<M; j++)
            set_B(i,j, PSEUDO);

        P[i]=PSEUDO;
        Q[i]=PSEUDO;
    }

    float64_t allpatprob=0 ;

#ifdef USE_HMMPARALLEL
    int32_t num_threads = parallel->get_num_threads();
    pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
    S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);

    if (p_observations->get_num_vectors()<num_threads)
        num_threads=p_observations->get_num_vectors();
#endif

    for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
    {

#ifdef USE_HMMPARALLEL
        if (dim%num_threads==0)
        {
            for (i=0; i<num_threads; i++)
            {
                if (dim+i<p_observations->get_num_vectors())
                {
                    params[i].hmm=estimate ;
                    params[i].dim=dim+i ;
                    pthread_create(&threads[i], NULL, vit_dim_prefetch, (void*)&params[i]) ;
                }
            }
            for (i=0; i<num_threads; i++)
            {
                if (dim+i<p_observations->get_num_vectors())
                {
                    pthread_join(threads[i], NULL);
                    allpatprob += params[i].prob_sum;
                }
            }
        }
#else
        //using viterbi to find best path
        allpatprob += estimate->best_path(dim);
#endif // USE_HMMPARALLEL

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

        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 );

        P[estimate->PATH(dim)[0]]++;
        Q[estimate->PATH(dim)[p_observations->get_vector_length(dim)-1]]++;
    }

#ifdef USE_HMMPARALLEL
    SG_FREE(threads);
    SG_FREE(params);
#endif

    allpatprob/=p_observations->get_num_vectors() ;
    estimate->all_pat_prob=allpatprob ;
    estimate->all_path_prob_updated=true ;

    //converting A to probability measure a
    for (i=0; i<N; i++)
    {
        sum=0;
        for (j=0; j<N; j++)
            sum+=get_A(i,j);

        for (j=0; j<N; j++)
            set_a(i,j, log(get_A(i,j)/sum));
    }

    //converting B to probability measures b
    for (i=0; i<N; i++)
    {
        sum=0;
        for (j=0; j<M; j++)
            sum+=get_B(i,j);

        for (j=0; j<M; j++)
            set_b(i,j, log(get_B(i, j)/sum));
    }

    //converting P to probability measure p
    sum=0;
    for (i=0; i<N; i++)
        sum+=P[i];

    for (i=0; i<N; i++)
        set_p(i, log(P[i]/sum));

    //converting Q to probability measure q
    sum=0;
    for (i=0; i<N; i++)
        sum+=Q[i];

    for (i=0; i<N; i++)
        set_q(i, log(Q[i]/sum));

    //new model probability is unknown
    invalidate_model();
}

// estimate parameters listed in learn_x
void CHMM::estimate_model_viterbi_defined(CHMM* estimate)
{
    int32_t i,j,k,t;
    float64_t sum;
    float64_t* P=ARRAYN1(0);
    float64_t* Q=ARRAYN2(0);

    path_deriv_updated=false ;

    //initialize with pseudocounts
    for (i=0; i<N; i++)
    {
        for (j=0; j<N; j++)
            set_A(i,j, PSEUDO);

        for (j=0; j<M; j++)
            set_B(i,j, PSEUDO);

        P[i]=PSEUDO;
        Q[i]=PSEUDO;
    }

#ifdef USE_HMMPARALLEL
    int32_t num_threads = parallel->get_num_threads();
    pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
    S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);
#endif

    float64_t allpatprob=0.0 ;
    for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
    {

#ifdef USE_HMMPARALLEL
        if (dim%num_threads==0)
        {
            for (i=0; i<num_threads; i++)
            {
                if (dim+i<p_observations->get_num_vectors())
                {
                    params[i].hmm=estimate ;
                    params[i].dim=dim+i ;
                    pthread_create(&threads[i], NULL, vit_dim_prefetch, (void*)&params[i]) ;
                }
            }
            for (i=0; i<num_threads; i++)
            {
                if (dim+i<p_observations->get_num_vectors())
                {
                    pthread_join(threads[i], NULL);
                    allpatprob += params[i].prob_sum;
                }
            }
        }
#else // USE_HMMPARALLEL
        //using viterbi to find best path
        allpatprob += estimate->best_path(dim);
#endif // USE_HMMPARALLEL


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

        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 );

        P[estimate->PATH(dim)[0]]++;
        Q[estimate->PATH(dim)[p_observations->get_vector_length(dim)-1]]++;
    }

#ifdef USE_HMMPARALLEL
    SG_FREE(threads);
    SG_FREE(params);
#endif

    //estimate->invalidate_model() ;
    //float64_t q=estimate->best_path(-1) ;

    allpatprob/=p_observations->get_num_vectors() ;
    estimate->all_pat_prob=allpatprob ;
    estimate->all_path_prob_updated=true ;


    //copy old model
    for (i=0; i<N; i++)
    {
        for (j=0; j<N; j++)
            set_a(i,j, estimate->get_a(i,j));

        for (j=0; j<M; j++)
            set_b(i,j, estimate->get_b(i,j));
    }

    //converting A to probability measure a
    i=0;
    sum=0;
    j=model->get_learn_a(i,0);
    k=i;
    while (model->get_learn_a(i,0)!=-1 || k<i)
    {
        if (j==model->get_learn_a(i,0))
        {
            sum+=get_A(model->get_learn_a(i,0), model->get_learn_a(i,1));
            i++;
        }
        else
        {
            while (k<i)
            {
                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));
                k++;
            }

            sum=0;
            j=model->get_learn_a(i,0);
            k=i;
        }
    }

    //converting B to probability measures b
    i=0;
    sum=0;
    j=model->get_learn_b(i,0);
    k=i;
    while (model->get_learn_b(i,0)!=-1 || k<i)
    {
        if (j==model->get_learn_b(i,0))
        {
            sum+=get_B(model->get_learn_b(i,0),model->get_learn_b(i,1));
            i++;
        }
        else
        {
            while (k<i)
            {
                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));
                k++;
            }

            sum=0;
            j=model->get_learn_b(i,0);
            k=i;
        }
    }

    i=0;
    sum=0;
    while (model->get_learn_p(i)!=-1)
    {
        sum+=P[model->get_learn_p(i)] ;
        i++ ;
    } ;
    i=0 ;
    while (model->get_learn_p(i)!=-1)
    {
        set_p(model->get_learn_p(i), log(P[model->get_learn_p(i)]/sum));
        i++ ;
    } ;

    i=0;
    sum=0;
    while (model->get_learn_q(i)!=-1)
    {
        sum+=Q[model->get_learn_q(i)] ;
        i++ ;
    } ;
    i=0 ;
    while (model->get_learn_q(i)!=-1)
    {
        set_q(model->get_learn_q(i), log(Q[model->get_learn_q(i)]/sum));
        i++ ;
    } ;


    //new model probability is unknown
    invalidate_model();
}
//------------------------------------------------------------------------------------//

//to give an idea what the model looks like
void CHMM::output_model(bool verbose)
{
    int32_t i,j;
    float64_t checksum;

    //generic info
    SG_INFO("log(Pr[O|model])=%e, #states: %i, #observationssymbols: %i, #observations: %ix%i\n",
            (float64_t)((p_observations) ? model_probability() : -CMath::INFTY),
            N, M, ((p_observations) ? p_observations->get_max_vector_length() : 0), ((p_observations) ? p_observations->get_num_vectors() : 0));

    if (verbose)
    {
        // tranisition matrix a
        SG_INFO("\ntransition matrix\n")
        for (i=0; i<N; i++)
        {
            checksum= get_q(i);
            for (j=0; j<N; j++)
            {
                checksum= CMath::logarithmic_sum(checksum, get_a(i,j));

                SG_INFO("a(%02i,%02i)=%1.4f ",i,j, (float32_t) exp(get_a(i,j)))

                if (j % 4 == 3)
                    SG_PRINT("\n")
            }
            if (fabs(checksum)>1e-5)
                SG_DEBUG(" checksum % E ******* \n",checksum)
            else
                SG_DEBUG(" checksum % E\n",checksum)
        }

        // distribution of start states p
        SG_INFO("\ndistribution of start states\n")
        checksum=-CMath::INFTY;
        for (i=0; i<N; i++)
        {
            checksum= CMath::logarithmic_sum(checksum, get_p(i));
            SG_INFO("p(%02i)=%1.4f ",i, (float32_t) exp(get_p(i)))
            if (i % 4 == 3)
                SG_PRINT("\n")
        }
        if (fabs(checksum)>1e-5)
            SG_DEBUG(" checksum % E ******* \n",checksum)
        else
            SG_DEBUG(" checksum=% E\n", checksum)

        // distribution of terminal states p
        SG_INFO("\ndistribution of terminal states\n")
        checksum=-CMath::INFTY;
        for (i=0; i<N; i++)
        {
            checksum= CMath::logarithmic_sum(checksum, get_q(i));
            SG_INFO("q(%02i)=%1.4f ",i, (float32_t) exp(get_q(i)))
            if (i % 4 == 3)
                SG_INFO("\n")
        }
        if (fabs(checksum)>1e-5)
            SG_DEBUG(" checksum % E ******* \n",checksum)
        else
            SG_DEBUG(" checksum=% E\n", checksum)

        // distribution of observations given the state b
        SG_INFO("\ndistribution of observations given the state\n")
        for (i=0; i<N; i++)
        {
            checksum=-CMath::INFTY;
            for (j=0; j<M; j++)
            {
                checksum=CMath::logarithmic_sum(checksum, get_b(i,j));
                SG_INFO("b(%02i,%02i)=%1.4f ",i,j, (float32_t) exp(get_b(i,j)))
                if (j % 4 == 3)
                    SG_PRINT("\n")
            }
            if (fabs(checksum)>1e-5)
                SG_DEBUG(" checksum % E ******* \n",checksum)
            else
                SG_DEBUG(" checksum % E\n",checksum)
        }
    }
    SG_PRINT("\n")
}

//to give an idea what the model looks like
void CHMM::output_model_defined(bool verbose)
{
    int32_t i,j;
    if (!model)
        return ;

    //generic info
    SG_INFO("log(Pr[O|model])=%e, #states: %i, #observationssymbols: %i, #observations: %ix%i\n",
            (float64_t)((p_observations) ? model_probability() : -CMath::INFTY),
            N, M, ((p_observations) ? p_observations->get_max_vector_length() : 0), ((p_observations) ? p_observations->get_num_vectors() : 0));

    if (verbose)
    {
        // tranisition matrix a
        SG_INFO("\ntransition matrix\n")

        //initialize a values that have to be learned
        i=0;
        j=model->get_learn_a(i,0);
        while (model->get_learn_a(i,0)!=-1)
        {
            if (j!=model->get_learn_a(i,0))
            {
                j=model->get_learn_a(i,0);
                SG_PRINT("\n")
            }

            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))))
            i++;
        }

        // distribution of observations given the state b
        SG_INFO("\n\ndistribution of observations given the state\n")
        i=0;
        j=model->get_learn_b(i,0);
        while (model->get_learn_b(i,0)!=-1)
        {
            if (j!=model->get_learn_b(i,0))
            {
                j=model->get_learn_b(i,0);
                SG_PRINT("\n")
            }

            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))))
            i++;
        }

        SG_PRINT("\n")
    }
    SG_PRINT("\n")
}

//------------------------------------------------------------------------------------//

//convert model to log probabilities
void CHMM::convert_to_log()
{
    int32_t i,j;

    for (i=0; i<N; i++)
    {
        if (get_p(i)!=0)
            set_p(i, log(get_p(i)));
        else
            set_p(i, -CMath::INFTY);;
    }

    for (i=0; i<N; i++)
    {
        if (get_q(i)!=0)
            set_q(i, log(get_q(i)));
        else
            set_q(i, -CMath::INFTY);;
    }

    for (i=0; i<N; i++)
    {
        for (j=0; j<N; j++)
        {
            if (get_a(i,j)!=0)
                set_a(i,j, log(get_a(i,j)));
            else
                set_a(i,j, -CMath::INFTY);
        }
    }

    for (i=0; i<N; i++)
    {
        for (j=0; j<M; j++)
        {
            if (get_b(i,j)!=0)
                set_b(i,j, log(get_b(i,j)));
            else
                set_b(i,j, -CMath::INFTY);
        }
    }
    loglikelihood=true;

    invalidate_model();
}

//init model with random values
void CHMM::init_model_random()
{
    const float64_t MIN_RAND=23e-3;

    float64_t sum;
    int32_t i,j;

    //initialize a with random values
    for (i=0; i<N; i++)
    {
        sum=0;
        for (j=0; j<N; j++)
        {
            set_a(i,j, CMath::random(MIN_RAND, 1.0));

            sum+=get_a(i,j);
        }

        for (j=0; j<N; j++)
            set_a(i,j, get_a(i,j)/sum);
    }

    //initialize pi with random values
    sum=0;
    for (i=0; i<N; i++)
    {
        set_p(i, CMath::random(MIN_RAND, 1.0));

        sum+=get_p(i);
    }

    for (i=0; i<N; i++)
        set_p(i, get_p(i)/sum);

    //initialize q with random values
    sum=0;
    for (i=0; i<N; i++)
    {
        set_q(i, CMath::random(MIN_RAND, 1.0));

        sum+=get_q(i);
    }

    for (i=0; i<N; i++)
        set_q(i, get_q(i)/sum);

    //initialize b with random values
    for (i=0; i<N; i++)
    {
        sum=0;
        for (j=0; j<M; j++)
        {
            set_b(i,j, CMath::random(MIN_RAND, 1.0));

            sum+=get_b(i,j);
        }

        for (j=0; j<M; j++)
            set_b(i,j, get_b(i,j)/sum);
    }

    //initialize pat/mod_prob as not calculated
    invalidate_model();
}

//init model according to const_x
void CHMM::init_model_defined()
{
    int32_t i,j,k,r;
    float64_t sum;
    const float64_t MIN_RAND=23e-3;

    //initialize a with zeros
    for (i=0; i<N; i++)
        for (j=0; j<N; j++)
            set_a(i,j, 0);

    //initialize p with zeros
    for (i=0; i<N; i++)
        set_p(i, 0);

    //initialize q with zeros
    for (i=0; i<N; i++)
        set_q(i, 0);

    //initialize b with zeros
    for (i=0; i<N; i++)
        for (j=0; j<M; j++)
            set_b(i,j, 0);


    //initialize a values that have to be learned
    float64_t *R=SG_MALLOC(float64_t, N);
    for (r=0; r<N; r++) R[r]=CMath::random(MIN_RAND,1.0);
    i=0; sum=0; k=i;
    j=model->get_learn_a(i,0);
    while (model->get_learn_a(i,0)!=-1 || k<i)
    {
        if (j==model->get_learn_a(i,0))
        {
            sum+=R[model->get_learn_a(i,1)] ;
            i++;
        }
        else
        {
            while (k<i)
            {
                set_a(model->get_learn_a(k,0), model->get_learn_a(k,1),
                        R[model->get_learn_a(k,1)]/sum);
                k++;
            }
            j=model->get_learn_a(i,0);
            k=i;
            sum=0;
            for (r=0; r<N; r++) R[r]=CMath::random(MIN_RAND,1.0);
        }
    }
    SG_FREE(R); R=NULL ;

    //initialize b values that have to be learned
    R=SG_MALLOC(float64_t, M);
    for (r=0; r<M; r++) R[r]=CMath::random(MIN_RAND,1.0);
    i=0; sum=0; k=0 ;
    j=model->get_learn_b(i,0);
    while (model->get_learn_b(i,0)!=-1 || k<i)
    {
        if (j==model->get_learn_b(i,0))
        {
            sum+=R[model->get_learn_b(i,1)] ;
            i++;
        }
        else
        {
            while (k<i)
            {
                set_b(model->get_learn_b(k,0),model->get_learn_b(k,1),
                        R[model->get_learn_b(k,1)]/sum);
                k++;
            }

            j=model->get_learn_b(i,0);
            k=i;
            sum=0;
            for (r=0; r<M; r++) R[r]=CMath::random(MIN_RAND,1.0);
        }
    }
    SG_FREE(R); R=NULL ;

    //set consts into a
    i=0;
    while (model->get_const_a(i,0) != -1)
    {
        set_a(model->get_const_a(i,0), model->get_const_a(i,1), model->get_const_a_val(i));
        i++;
    }

    //set consts into b
    i=0;
    while (model->get_const_b(i,0) != -1)
    {
        set_b(model->get_const_b(i,0), model->get_const_b(i,1), model->get_const_b_val(i));
        i++;
    }

    //set consts into p
    i=0;
    while (model->get_const_p(i) != -1)
    {
        set_p(model->get_const_p(i), model->get_const_p_val(i));
        i++;
    }

    //initialize q with zeros
    for (i=0; i<N; i++)
        set_q(i, 0.0);

    //set consts into q
    i=0;
    while (model->get_const_q(i) != -1)
    {
        set_q(model->get_const_q(i), model->get_const_q_val(i));
        i++;
    }

    // init p
    i=0;
    sum=0;
    while (model->get_learn_p(i)!=-1)
    {
        set_p(model->get_learn_p(i),CMath::random(MIN_RAND,1.0)) ;
        sum+=get_p(model->get_learn_p(i)) ;
        i++ ;
    } ;
    i=0 ;
    while (model->get_learn_p(i)!=-1)
    {
        set_p(model->get_learn_p(i), get_p(model->get_learn_p(i))/sum);
        i++ ;
    } ;

    // initialize q
    i=0;
    sum=0;
    while (model->get_learn_q(i)!=-1)
    {
        set_q(model->get_learn_q(i),CMath::random(MIN_RAND,1.0)) ;
        sum+=get_q(model->get_learn_q(i)) ;
        i++ ;
    } ;
    i=0 ;
    while (model->get_learn_q(i)!=-1)
    {
        set_q(model->get_learn_q(i), get_q(model->get_learn_q(i))/sum);
        i++ ;
    } ;

    //initialize pat/mod_prob as not calculated
    invalidate_model();
}

void CHMM::clear_model()
{
    int32_t i,j;
    for (i=0; i<N; i++)
    {
        set_p(i, log(PSEUDO));
        set_q(i, log(PSEUDO));

        for (j=0; j<N; j++)
            set_a(i,j, log(PSEUDO));

        for (j=0; j<M; j++)
            set_b(i,j, log(PSEUDO));
    }
}

void CHMM::clear_model_defined()
{
    int32_t i,j,k;

    for (i=0; (j=model->get_learn_p(i))!=-1; i++)
        set_p(j, log(PSEUDO));

    for (i=0; (j=model->get_learn_q(i))!=-1; i++)
        set_q(j, log(PSEUDO));

    for (i=0; (j=model->get_learn_a(i,0))!=-1; i++)
    {
        k=model->get_learn_a(i,1); // catch (j,k) as indizes to be learned
        set_a(j,k, log(PSEUDO));
    }

    for (i=0; (j=model->get_learn_b(i,0))!=-1; i++)
    {
        k=model->get_learn_b(i,1); // catch (j,k) as indizes to be learned
        set_b(j,k, log(PSEUDO));
    }
}

void CHMM::copy_model(CHMM* l)
{
    int32_t i,j;
    for (i=0; i<N; i++)
    {
        set_p(i, l->get_p(i));
        set_q(i, l->get_q(i));

        for (j=0; j<N; j++)
            set_a(i,j, l->get_a(i,j));

        for (j=0; j<M; j++)
            set_b(i,j, l->get_b(i,j));
    }
}

void CHMM::invalidate_model()
{
    //initialize pat/mod_prob/alpha/beta cache as not calculated
    this->mod_prob=0.0;
    this->mod_prob_updated=false;

    if (mem_initialized)
    {
      if (trans_list_forward_cnt)
        SG_FREE(trans_list_forward_cnt);
      trans_list_forward_cnt=NULL ;
      if (trans_list_backward_cnt)
        SG_FREE(trans_list_backward_cnt);
      trans_list_backward_cnt=NULL ;
      if (trans_list_forward)
        {
          for (int32_t i=0; i<trans_list_len; i++)
        if (trans_list_forward[i])
          SG_FREE(trans_list_forward[i]);
          SG_FREE(trans_list_forward);
          trans_list_forward=NULL ;
        }
      if (trans_list_backward)
        {
          for (int32_t i=0; i<trans_list_len; i++)
        if (trans_list_backward[i])
          SG_FREE(trans_list_backward[i]);
          SG_FREE(trans_list_backward);
          trans_list_backward = NULL ;
        } ;

      trans_list_len = N ;
      trans_list_forward = SG_MALLOC(T_STATES*, N);
      trans_list_forward_cnt = SG_MALLOC(T_STATES, N);

      for (int32_t j=0; j<N; j++)
        {
          trans_list_forward_cnt[j]= 0 ;
          trans_list_forward[j]= SG_MALLOC(T_STATES, N);
          for (int32_t i=0; i<N; i++)
        if (get_a(i,j)>CMath::ALMOST_NEG_INFTY)
          {
            trans_list_forward[j][trans_list_forward_cnt[j]]=i ;
            trans_list_forward_cnt[j]++ ;
          }
        } ;

      trans_list_backward = SG_MALLOC(T_STATES*, N);
      trans_list_backward_cnt = SG_MALLOC(T_STATES, N);

      for (int32_t i=0; i<N; i++)
        {
          trans_list_backward_cnt[i]= 0 ;
          trans_list_backward[i]= SG_MALLOC(T_STATES, N);
          for (int32_t j=0; j<N; j++)
        if (get_a(i,j)>CMath::ALMOST_NEG_INFTY)
          {
            trans_list_backward[i][trans_list_backward_cnt[i]]=j ;
            trans_list_backward_cnt[i]++ ;
          }
        } ;
    } ;
    this->all_pat_prob=0.0;
    this->pat_prob=0.0;
    this->path_deriv_updated=false ;
    this->path_deriv_dimension=-1 ;
    this->all_path_prob_updated=false;

#ifdef USE_HMMPARALLEL_STRUCTURES
    {
        for (int32_t i=0; i<parallel->get_num_threads(); i++)
        {
            this->alpha_cache[i].updated=false;
            this->beta_cache[i].updated=false;
            path_prob_updated[i]=false ;
            path_prob_dimension[i]=-1 ;
        } ;
    }
#else // USE_HMMPARALLEL_STRUCTURES
    this->alpha_cache.updated=false;
    this->beta_cache.updated=false;
    this->path_prob_dimension=-1;
    this->path_prob_updated=false;

#endif // USE_HMMPARALLEL_STRUCTURES
}

void CHMM::open_bracket(FILE* file)
{
    int32_t value;
    while (((value=fgetc(file)) != EOF) && (value!='['))    //skip possible spaces and end if '[' occurs
    {
        if (value=='\n')
            line++;
    }

    if (value==EOF)
        error(line, "expected \"[\" in input file");

    while (((value=fgetc(file)) != EOF) && (isspace(value)))    //skip possible spaces
    {
        if (value=='\n')
            line++;
    }

    ungetc(value, file);
}

void CHMM::close_bracket(FILE* file)
{
    int32_t value;
    while (((value=fgetc(file)) != EOF) && (value!=']'))    //skip possible spaces and end if ']' occurs
    {
        if (value=='\n')
            line++;
    }

    if (value==EOF)
        error(line, "expected \"]\" in input file");
}

bool CHMM::comma_or_space(FILE* file)
{
    int32_t value;
    while (((value=fgetc(file)) != EOF) && (value!=',') && (value!=';') && (value!=']'))     //skip possible spaces and end if ',' or ';' occurs
    {
        if (value=='\n')
            line++;
    }
    if (value==']')
    {
        ungetc(value, file);
        SG_ERROR("found ']' instead of ';' or ','\n")
        return false ;
    } ;

    if (value==EOF)
        error(line, "expected \";\" or \",\" in input file");

    while (((value=fgetc(file)) != EOF) && (isspace(value)))    //skip possible spaces
    {
        if (value=='\n')
            line++;
    }
    ungetc(value, file);
    return true ;
}

bool CHMM::get_numbuffer(FILE* file, char* buffer, int32_t length)
{
    signed char value;

    while (((value=fgetc(file)) != EOF) &&
            !isdigit(value) && (value!='A')
            && (value!='C') && (value!='G') && (value!='T')
            && (value!='N') && (value!='n')
            && (value!='.') && (value!='-') && (value!='e') && (value!=']')) //skip possible spaces+crap
    {
        if (value=='\n')
            line++;
    }
    if (value==']')
    {
        ungetc(value,file) ;
        return false ;
    } ;
    if (value!=EOF)
    {
        int32_t i=0;
        switch (value)
        {
            case 'A':
                value='0' +CAlphabet::B_A;
                break;
            case 'C':
                value='0' +CAlphabet::B_C;
                break;
            case 'G':
                value='0' +CAlphabet::B_G;
                break;
            case 'T':
                value='0' +CAlphabet::B_T;
                break;
        };

        buffer[i++]=value;

        while (((value=fgetc(file)) != EOF) &&
                (isdigit(value) || (value=='.') || (value=='-') || (value=='e')
                 || (value=='A') || (value=='C') || (value=='G')|| (value=='T')
                 || (value=='N') || (value=='n')) && (i<length))
        {
            switch (value)
            {
                case 'A':
                    value='0' +CAlphabet::B_A;
                    break;
                case 'C':
                    value='0' +CAlphabet::B_C;
                    break;
                case 'G':
                    value='0' +CAlphabet::B_G;
                    break;
                case 'T':
                    value='0' +CAlphabet::B_T;
                    break;
                case '1': case '2': case'3': case '4': case'5':
                case '6': case '7': case'8': case '9': case '0': break ;
                case '.': case 'e': case '-': break ;
                default:
                                              SG_ERROR("found crap: %i %c (pos:%li)\n",i,value,ftell(file))
            };
            buffer[i++]=value;
        }
        ungetc(value, file);
        buffer[i]='\0';

        return (i<=length) && (i>0);
    }
    return false;
}

/*
   -format specs: model_file (model.hmm)
   % HMM - specification
   % N  - number of states
   % M  - number of observation_tokens
   % a is state_transition_matrix
   % size(a)= [N,N]
   %
   % b is observation_per_state_matrix
   % size(b)= [N,M]
   %
   % p is initial distribution
   % size(p)= [1, N]

   N=<int32_t>;
   M=<int32_t>;

   p=[<float64_t>,<float64_t>...<DOUBLE>];
   q=[<DOUBLE>,<DOUBLE>...<DOUBLE>];

   a=[ [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   ];

   b=[ [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   ];
   */

bool CHMM::load_model(FILE* file)
{
    int32_t received_params=0;  //a,b,p,N,M,O

    bool result=false;
    E_STATE state=INITIAL;
    char buffer[1024];

    line=1;
    int32_t i,j;

    if (file)
    {
        while (state!=END)
        {
            int32_t value=fgetc(file);

            if (value=='\n')
                line++;
            if (value==EOF)
                state=END;

            switch (state)
            {
                case INITIAL:   // in the initial state only N,M initialisations and comments are allowed
                    if (value=='N')
                    {
                        if (received_params & GOTN)
                            error(line, "in model file: \"p double defined\"");
                        else
                            state=GET_N;
                    }
                    else if (value=='M')
                    {
                        if (received_params & GOTM)
                            error(line, "in model file: \"p double defined\"");
                        else
                            state=GET_M;
                    }
                    else if (value=='%')
                    {
                        state=COMMENT;
                    }
                    break;
                case ARRAYs:    // when n,m, order are known p,a,b arrays are allowed to be read
                    if (value=='p')
                    {
                        if (received_params & GOTp)
                            error(line, "in model file: \"p double defined\"");
                        else
                            state=GET_p;
                    }
                    if (value=='q')
                    {
                        if (received_params & GOTq)
                            error(line, "in model file: \"q double defined\"");
                        else
                            state=GET_q;
                    }
                    else if (value=='a')
                    {
                        if (received_params & GOTa)
                            error(line, "in model file: \"a double defined\"");
                        else
                            state=GET_a;
                    }
                    else if (value=='b')
                    {
                        if (received_params & GOTb)
                            error(line, "in model file: \"b double defined\"");
                        else
                            state=GET_b;
                    }
                    else if (value=='%')
                    {
                        state=COMMENT;
                    }
                    break;
                case GET_N:
                    if (value=='=')
                    {
                        if (get_numbuffer(file, buffer, 4)) //get num
                        {
                            this->N= atoi(buffer);
                            received_params|=GOTN;
                            state= (received_params == (GOTN | GOTM | GOTO)) ? ARRAYs : INITIAL;
                        }
                        else
                            state=END;      //end if error
                    }
                    break;
                case GET_M:
                    if (value=='=')
                    {
                        if (get_numbuffer(file, buffer, 4)) //get num
                        {
                            this->M= atoi(buffer);
                            received_params|=GOTM;
                            state= (received_params == (GOTN | GOTM | GOTO)) ? ARRAYs : INITIAL;
                        }
                        else
                            state=END;      //end if error
                    }
                    break;
                case GET_a:
                    if (value=='=')
                    {
                        float64_t f;

                        transition_matrix_a=SG_MALLOC(float64_t, N*N);
                        open_bracket(file);
                        for (i=0; i<this->N; i++)
                        {
                            open_bracket(file);

                            for (j=0; j<this->N ; j++)
                            {

                                if (fscanf( file, "%le", &f ) != 1)
                                    error(line, "float64_t expected");
                                else
                                    set_a(i,j, f);

                                if (j<this->N-1)
                                    comma_or_space(file);
                                else
                                    close_bracket(file);
                            }

                            if (i<this->N-1)
                                comma_or_space(file);
                            else
                                close_bracket(file);
                        }
                        received_params|=GOTa;
                    }
                    state= (received_params == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
                    break;
                case GET_b:
                    if (value=='=')
                    {
                        float64_t f;

                        observation_matrix_b=SG_MALLOC(float64_t, N*M);
                        open_bracket(file);
                        for (i=0; i<this->N; i++)
                        {
                            open_bracket(file);

                            for (j=0; j<this->M ; j++)
                            {

                                if (fscanf( file, "%le", &f ) != 1)
                                    error(line, "float64_t expected");
                                else
                                    set_b(i,j, f);

                                if (j<this->M-1)
                                    comma_or_space(file);
                                else
                                    close_bracket(file);
                            }

                            if (i<this->N-1)
                                comma_or_space(file);
                            else
                                close_bracket(file);
                        }
                        received_params|=GOTb;
                    }
                    state= ((received_params & (GOTa | GOTb | GOTp | GOTq)) == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
                    break;
                case GET_p:
                    if (value=='=')
                    {
                        float64_t f;

                        initial_state_distribution_p=SG_MALLOC(float64_t, N);
                        open_bracket(file);
                        for (i=0; i<this->N ; i++)
                        {
                            if (fscanf( file, "%le", &f ) != 1)
                                error(line, "float64_t expected");
                            else
                                set_p(i, f);

                            if (i<this->N-1)
                                comma_or_space(file);
                            else
                                close_bracket(file);
                        }
                        received_params|=GOTp;
                    }
                    state= (received_params == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
                    break;
                case GET_q:
                    if (value=='=')
                    {
                        float64_t f;

                        end_state_distribution_q=SG_MALLOC(float64_t, N);
                        open_bracket(file);
                        for (i=0; i<this->N ; i++)
                        {
                            if (fscanf( file, "%le", &f ) != 1)
                                error(line, "float64_t expected");
                            else
                                set_q(i, f);

                            if (i<this->N-1)
                                comma_or_space(file);
                            else
                                close_bracket(file);
                        }
                        received_params|=GOTq;
                    }
                    state= (received_params == (GOTa | GOTb | GOTp | GOTq)) ? END : ARRAYs;
                    break;
                case COMMENT:
                    if (value==EOF)
                        state=END;
                    else if (value=='\n')
                    {
                        line++;
                        state=INITIAL;
                    }
                    break;

                default:
                    break;
            }
        }
        result= (received_params== (GOTa | GOTb | GOTp | GOTq | GOTN | GOTM | GOTO));
    }

    SG_WARNING("not normalizing anymore, call normalize_hmm to make sure the hmm is valid!!\n")
    return result;
}

/*
    -format specs: train_file (train.trn)
    % HMM-TRAIN - specification
    % learn_a - elements in state_transition_matrix to be learned
    % learn_b - elements in oberservation_per_state_matrix to be learned
    %           note: each line stands for
    %               <state>, <observation(0)>, observation(1)...observation(NOW)>
    % learn_p - elements in initial distribution to be learned
    % learn_q - elements in the end-state distribution to be learned
    %
    % const_x - specifies initial values of elements
    %               rest is assumed to be 0.0
    %
    %   NOTE: IMPLICIT DEFINES:
    %       #define A 0
    %       #define C 1
    %       #define G 2
    %       #define T 3
    %

    learn_a=[ [<int32_t>,<int32_t>];
    [<int32_t>,<int32_t>];
    [<int32_t>,<int32_t>];
    ........
    [<int32_t>,<int32_t>];
    [-1,-1];
    ];

    learn_b=[ [<int32_t>,<int32_t>];
    [<int32_t>,<int32_t>];
    [<int32_t>,<int32_t>];
    ........
    [<int32_t>,<int32_t>];
    [-1,-1];
    ];

    learn_p= [ <int32_t>, ... , <int32_t>, -1 ];
    learn_q= [ <int32_t>, ... , <int32_t>, -1 ];


    const_a=[ [<int32_t>,<int32_t>,<DOUBLE>];
    [<int32_t>,<int32_t>,<DOUBLE>];
    [<int32_t>,<int32_t>,<DOUBLE>];
    ........
    [<int32_t>,<int32_t>,<DOUBLE>];
    [-1,-1,-1];
    ];

    const_b=[ [<int32_t>,<int32_t>,<DOUBLE>];
    [<int32_t>,<int32_t>,<DOUBLE>];
    [<int32_t>,<int32_t>,<DOUBLE];
    ........
    [<int32_t>,<int32_t>,<DOUBLE>];
    [-1,-1];
    ];

    const_p[]=[ [<int32_t>, <DOUBLE>], ... , [<int32_t>,<DOUBLE>], [-1,-1] ];
    const_q[]=[ [<int32_t>, <DOUBLE>], ... , [<int32_t>,<DOUBLE>], [-1,-1] ];
    */
bool CHMM::load_definitions(FILE* file, bool verbose, bool _initialize)
{
    if (model)
        delete model ;
    model=new Model();

    int32_t received_params=0x0000000;  //a,b,p,q,N,M,O
    char buffer[1024];

    bool result=false;
    E_STATE state=INITIAL;

    { // do some useful initializations
        model->set_learn_a(0, -1);
        model->set_learn_a(1, -1);
        model->set_const_a(0, -1);
        model->set_const_a(1, -1);
        model->set_const_a_val(0, 1.0);
        model->set_learn_b(0, -1);
        model->set_const_b(0, -1);
        model->set_const_b_val(0, 1.0);
        model->set_learn_p(0, -1);
        model->set_learn_q(0, -1);
        model->set_const_p(0, -1);
        model->set_const_q(0, -1);
    } ;

    line=1;

    if (file)
    {
        while (state!=END)
        {
            int32_t value=fgetc(file);

            if (value=='\n')
                line++;

            if (value==EOF)
                state=END;

            switch (state)
            {
                case INITIAL:
                    if (value=='l')
                    {
                        if (fgetc(file)=='e' && fgetc(file)=='a' && fgetc(file)=='r' && fgetc(file)=='n' && fgetc(file)=='_')
                        {
                            switch(fgetc(file))
                            {
                                case 'a':
                                    state=GET_learn_a;
                                    break;
                                case 'b':
                                    state=GET_learn_b;
                                    break;
                                case 'p':
                                    state=GET_learn_p;
                                    break;
                                case 'q':
                                    state=GET_learn_q;
                                    break;
                                default:
                                    error(line, "a,b,p or q expected in train definition file");
                            };
                        }
                    }
                    else if (value=='c')
                    {
                        if (fgetc(file)=='o' && fgetc(file)=='n' && fgetc(file)=='s'
                                && fgetc(file)=='t' && fgetc(file)=='_')
                        {
                            switch(fgetc(file))
                            {
                                case 'a':
                                    state=GET_const_a;
                                    break;
                                case 'b':
                                    state=GET_const_b;
                                    break;
                                case 'p':
                                    state=GET_const_p;
                                    break;
                                case 'q':
                                    state=GET_const_q;
                                    break;
                                default:
                                    error(line, "a,b,p or q expected in train definition file");
                            };
                        }
                    }
                    else if (value=='%')
                    {
                        state=COMMENT;
                    }
                    else if (value==EOF)
                    {
                        state=END;
                    }
                    break;
                case GET_learn_a:
                    if (value=='=')
                    {
                        open_bracket(file);
                        bool finished=false;
                        int32_t i=0;

                        if (verbose)
                            SG_DEBUG("\nlearn for transition matrix: ")
                        while (!finished)
                        {
                            open_bracket(file);

                            if (get_numbuffer(file, buffer, 4)) //get num
                            {
                                value=atoi(buffer);
                                model->set_learn_a(i++, value);

                                if (value<0)
                                {
                                    finished=true;
                                    break;
                                }
                                if (value>=N)
                                    SG_ERROR("invalid value for learn_a(%i,0): %i\n",i/2,(int)value)
                            }
                            else
                                break;

                            comma_or_space(file);

                            if (get_numbuffer(file, buffer, 4)) //get num
                            {
                                value=atoi(buffer);
                                model->set_learn_a(i++, value);

                                if (value<0)
                                {
                                    finished=true;
                                    break;
                                }
                                if (value>=N)
                                    SG_ERROR("invalid value for learn_a(%i,1): %i\n",i/2-1,(int)value)

                            }
                            else
                                break;
                            close_bracket(file);
                        }
                        close_bracket(file);
                        if (verbose)
                            SG_DEBUG("%i Entries",(int)(i/2))

                        if (finished)
                        {
                            received_params|=GOTlearn_a;

                            state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
                        }
                        else
                            state=END;
                    }
                    break;
                case GET_learn_b:
                    if (value=='=')
                    {
                        open_bracket(file);
                        bool finished=false;
                        int32_t i=0;

                        if (verbose)
                            SG_DEBUG("\nlearn for emission matrix:   ")

                        while (!finished)
                        {
                            open_bracket(file);

                            int32_t combine=0;

                            for (int32_t j=0; j<2; j++)
                            {
                                if (get_numbuffer(file, buffer, 4))   //get num
                                {
                                    value=atoi(buffer);

                                    if (j==0)
                                    {
                                        model->set_learn_b(i++, value);

                                        if (value<0)
                                        {
                                            finished=true;
                                            break;
                                        }
                                        if (value>=N)
                                            SG_ERROR("invalid value for learn_b(%i,0): %i\n",i/2,(int)value)
                                    }
                                    else
                                        combine=value;
                                }
                                else
                                    break;

                                if (j<1)
                                    comma_or_space(file);
                                else
                                    close_bracket(file);
                            }
                            model->set_learn_b(i++, combine);
                            if (combine>=M)

                                SG_ERROR("invalid value for learn_b(%i,1): %i\n",i/2-1,(int)value)
                        }
                        close_bracket(file);
                        if (verbose)
                            SG_DEBUG("%i Entries",(int)(i/2-1))

                        if (finished)
                        {
                            received_params|=GOTlearn_b;
                            state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
                        }
                        else
                            state=END;
                    }
                    break;
                case GET_learn_p:
                    if (value=='=')
                    {
                        open_bracket(file);
                        bool finished=false;
                        int32_t i=0;

                        if (verbose)
                            SG_DEBUG("\nlearn start states: ")
                        while (!finished)
                        {
                            if (get_numbuffer(file, buffer, 4)) //get num
                            {
                                value=atoi(buffer);

                                model->set_learn_p(i++, value);

                                if (value<0)
                                {
                                    finished=true;
                                    break;
                                }
                                if (value>=N)
                                    SG_ERROR("invalid value for learn_p(%i): %i\n",i-1,(int)value)
                            }
                            else
                                break;

                            comma_or_space(file);
                        }

                        close_bracket(file);
                        if (verbose)
                            SG_DEBUG("%i Entries",i-1)

                        if (finished)
                        {
                            received_params|=GOTlearn_p;
                            state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
                        }
                        else
                            state=END;
                    }
                    break;
                case GET_learn_q:
                    if (value=='=')
                    {
                        open_bracket(file);
                        bool finished=false;
                        int32_t i=0;

                        if (verbose)
                            SG_DEBUG("\nlearn terminal states: ")
                        while (!finished)
                        {
                            if (get_numbuffer(file, buffer, 4)) //get num
                            {
                                value=atoi(buffer);
                                model->set_learn_q(i++, value);

                                if (value<0)
                                {
                                    finished=true;
                                    break;
                                }
                                if (value>=N)
                                    SG_ERROR("invalid value for learn_q(%i): %i\n",i-1,(int)value)
                            }
                            else
                                break;

                            comma_or_space(file);
                        }

                        close_bracket(file);
                        if (verbose)
                            SG_DEBUG("%i Entries",i-1)

                        if (finished)
                        {
                            received_params|=GOTlearn_q;
                            state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
                        }
                        else
                            state=END;
                    }
                    break;
                case GET_const_a:
                    if (value=='=')
                    {
                        open_bracket(file);
                        bool finished=false;
                        int32_t i=0;

                        if (verbose)
#ifdef USE_HMMDEBUG
                            SG_DEBUG("\nconst for transition matrix: \n")
#else
                        SG_DEBUG("\nconst for transition matrix: ")
#endif
                        while (!finished)
                        {
                            open_bracket(file);

                            if (get_numbuffer(file, buffer, 4)) //get num
                            {
                                value=atoi(buffer);
                                model->set_const_a(i++, value);

                                if (value<0)
                                {
                                    finished=true;
                                    model->set_const_a(i++, value);
                                    model->set_const_a_val((int32_t)i/2 - 1, value);
                                    break;
                                }
                                if (value>=N)
                                    SG_ERROR("invalid value for const_a(%i,0): %i\n",i/2,(int)value)
                            }
                            else
                                break;

                            comma_or_space(file);

                            if (get_numbuffer(file, buffer, 4)) //get num
                            {
                                value=atoi(buffer);
                                model->set_const_a(i++, value);

                                if (value<0)
                                {
                                    finished=true;
                                    model->set_const_a_val((int32_t)i/2 - 1, value);
                                    break;
                                }
                                if (value>=N)
                                    SG_ERROR("invalid value for const_a(%i,1): %i\n",i/2-1,(int)value)
                            }
                            else
                                break;

                            if (!comma_or_space(file))
                                model->set_const_a_val((int32_t)i/2 - 1, 1.0);
                            else
                                if (get_numbuffer(file, buffer, 10))    //get num
                                {
                                    float64_t dvalue=atof(buffer);
                                    model->set_const_a_val((int32_t)i/2 - 1, dvalue);
                                    if (dvalue<0)
                                    {
                                        finished=true;
                                        break;
                                    }
                                    if ((dvalue>1.0) || (dvalue<0.0))
                                        SG_ERROR("invalid value for const_a_val(%i): %e\n",(int)i/2-1,dvalue)
                                }
                                else
                                    model->set_const_a_val((int32_t)i/2 - 1, 1.0);

#ifdef USE_HMMDEBUG
                            if (verbose)
                                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))
#endif
                            close_bracket(file);
                        }
                        close_bracket(file);
                        if (verbose)
                            SG_DEBUG("%i Entries",(int)i/2-1)

                        if (finished)
                        {
                            received_params|=GOTconst_a;
                            state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
                        }
                        else
                            state=END;
                    }
                    break;

                case GET_const_b:
                    if (value=='=')
                    {
                        open_bracket(file);
                        bool finished=false;
                        int32_t i=0;

                        if (verbose)
#ifdef USE_HMMDEBUG
                            SG_DEBUG("\nconst for emission matrix:   \n")
#else
                        SG_DEBUG("\nconst for emission matrix:   ")
#endif
                        while (!finished)
                        {
                            open_bracket(file);
                            int32_t combine=0;
                            for (int32_t j=0; j<3; j++)
                            {
                                if (get_numbuffer(file, buffer, 10))    //get num
                                {
                                    if (j==0)
                                    {
                                        value=atoi(buffer);

                                        model->set_const_b(i++, value);

                                        if (value<0)
                                        {
                                            finished=true;
                                            //model->set_const_b_val((int32_t)(i-1)/2, value);
                                            break;
                                        }
                                        if (value>=N)
                                            SG_ERROR("invalid value for const_b(%i,0): %i\n",i/2-1,(int)value)
                                    }
                                    else if (j==2)
                                    {
                                        float64_t dvalue=atof(buffer);
                                        model->set_const_b_val((int32_t)(i-1)/2, dvalue);
                                        if (dvalue<0)
                                        {
                                            finished=true;
                                            break;
                                        } ;
                                        if ((dvalue>1.0) || (dvalue<0.0))
                                            SG_ERROR("invalid value for const_b_val(%i,1): %e\n",i/2-1,dvalue)
                                    }
                                    else
                                    {
                                        value=atoi(buffer);
                                        combine= value;
                                    } ;
                                }
                                else
                                {
                                    if (j==2)
                                        model->set_const_b_val((int32_t)(i-1)/2, 1.0);
                                    break;
                                } ;
                                if (j<2)
                                    if ((!comma_or_space(file)) && (j==1))
                                    {
                                        model->set_const_b_val((int32_t)(i-1)/2, 1.0) ;
                                        break ;
                                    } ;
                            }
                            close_bracket(file);
                            model->set_const_b(i++, combine);
                            if (combine>=M)
                                SG_ERROR("invalid value for const_b(%i,1): %i\n",i/2-1, combine)
#ifdef USE_HMMDEBUG
                            if (verbose && !finished)
                                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))
#endif
                        }
                        close_bracket(file);
                        if (verbose)
                            SG_ERROR("%i Entries",(int)i/2-1)

                        if (finished)
                        {
                            received_params|=GOTconst_b;
                            state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
                        }
                        else
                            state=END;
                    }
                    break;
                case GET_const_p:
                    if (value=='=')
                    {
                        open_bracket(file);
                        bool finished=false;
                        int32_t i=0;

                        if (verbose)
#ifdef USE_HMMDEBUG
                            SG_DEBUG("\nconst for start states:     \n")
#else
                        SG_DEBUG("\nconst for start states:     ")
#endif
                        while (!finished)
                        {
                            open_bracket(file);

                            if (get_numbuffer(file, buffer, 4)) //get num
                            {
                                value=atoi(buffer);
                                model->set_const_p(i, value);

                                if (value<0)
                                {
                                    finished=true;
                                    model->set_const_p_val(i++, value);
                                    break;
                                }
                                if (value>=N)
                                    SG_ERROR("invalid value for const_p(%i): %i\n",i,(int)value)

                            }
                            else
                                break;

                            if (!comma_or_space(file))
                                model->set_const_p_val(i++, 1.0);
                            else
                                if (get_numbuffer(file, buffer, 10))    //get num
                                {
                                    float64_t dvalue=atof(buffer);
                                    model->set_const_p_val(i++, dvalue);
                                    if (dvalue<0)
                                    {
                                        finished=true;
                                        break;
                                    }
                                    if ((dvalue>1) || (dvalue<0))
                                        SG_ERROR("invalid value for const_p_val(%i): %e\n",i,dvalue)
                                }
                                else
                                    model->set_const_p_val(i++, 1.0);

                            close_bracket(file);

#ifdef USE_HMMDEBUG
                            if (verbose)
                                SG_DEBUG("const_p(%i)=%e\n", model->get_const_p(i-1),model->get_const_p_val(i-1))
#endif
                        }
                        if (verbose)
                            SG_DEBUG("%i Entries",i-1)

                        close_bracket(file);

                        if (finished)
                        {
                            received_params|=GOTconst_p;
                            state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
                        }
                        else
                            state=END;
                    }
                    break;
                case GET_const_q:
                    if (value=='=')
                    {
                        open_bracket(file);
                        bool finished=false;
                        if (verbose)
#ifdef USE_HMMDEBUG
                            SG_DEBUG("\nconst for terminal states: \n")
#else
                        SG_DEBUG("\nconst for terminal states: ")
#endif
                        int32_t i=0;

                        while (!finished)
                        {
                            open_bracket(file) ;
                            if (get_numbuffer(file, buffer, 4)) //get num
                            {
                                value=atoi(buffer);
                                model->set_const_q(i, value);
                                if (value<0)
                                {
                                    finished=true;
                                    model->set_const_q_val(i++, value);
                                    break;
                                }
                                if (value>=N)
                                    SG_ERROR("invalid value for const_q(%i): %i\n",i,(int)value)
                            }
                            else
                                break;

                            if (!comma_or_space(file))
                                model->set_const_q_val(i++, 1.0);
                            else
                                if (get_numbuffer(file, buffer, 10))    //get num
                                {
                                    float64_t dvalue=atof(buffer);
                                    model->set_const_q_val(i++, dvalue);
                                    if (dvalue<0)
                                    {
                                        finished=true;
                                        break;
                                    }
                                    if ((dvalue>1) || (dvalue<0))
                                        SG_ERROR("invalid value for const_q_val(%i): %e\n",i,(double) dvalue)
                                }
                                else
                                    model->set_const_q_val(i++, 1.0);

                            close_bracket(file);
#ifdef USE_HMMDEBUG
                            if (verbose)
                                SG_DEBUG("const_q(%i)=%e\n", model->get_const_q(i-1),model->get_const_q_val(i-1))
#endif
                        }
                        if (verbose)
                            SG_DEBUG("%i Entries",i-1)

                        close_bracket(file);

                        if (finished)
                        {
                            received_params|=GOTconst_q;
                            state= (received_params == (GOTlearn_a | GOTlearn_b | GOTlearn_p | GOTlearn_q |GOTconst_a | GOTconst_b | GOTconst_p | GOTconst_q)) ? END : INITIAL;
                        }
                        else
                            state=END;
                    }
                    break;
                case COMMENT:
                    if (value==EOF)
                        state=END;
                    else if (value=='\n')
                        state=INITIAL;
                    break;

                default:
                    break;
            }
        }
    }

    /*result=((received_params&(GOTlearn_a | GOTconst_a))!=0) ;
      result=result && ((received_params&(GOTlearn_b | GOTconst_b))!=0) ;
      result=result && ((received_params&(GOTlearn_p | GOTconst_p))!=0) ;
      result=result && ((received_params&(GOTlearn_q | GOTconst_q))!=0) ; */
    result=1 ;
    if (result)
    {
        model->sort_learn_a() ;
        model->sort_learn_b() ;
        if (_initialize)
        {
            init_model_defined(); ;
            convert_to_log();
        } ;
    }
    if (verbose)
        SG_DEBUG("\n")
    return result;
}

/*
   -format specs: model_file (model.hmm)
   % HMM - specification
   % N  - number of states
   % M  - number of observation_tokens
   % a is state_transition_matrix
   % size(a)= [N,N]
   %
   % b is observation_per_state_matrix
   % size(b)= [N,M]
   %
   % p is initial distribution
   % size(p)= [1, N]

   N=<int32_t>;
   M=<int32_t>;

   p=[<DOUBLE>,<DOUBLE>...<DOUBLE>];

   a=[ [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   ];

   b=[ [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   [<DOUBLE>,<DOUBLE>...<DOUBLE>];
   ];
   */

bool CHMM::save_model(FILE* file)
{
    bool result=false;
    int32_t i,j;
    const float32_t NAN_REPLACEMENT = (float32_t) CMath::ALMOST_NEG_INFTY ;

    if (file)
    {
        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");
        fprintf(file,"N=%d;\n",N);
        fprintf(file,"M=%d;\n",M);

        fprintf(file,"p=[");
        for (i=0; i<N; i++)
        {
            if (i<N-1) {
                if (CMath::is_finite(get_p(i)))
                    fprintf(file, "%e,", (double)get_p(i));
                else
                    fprintf(file, "%f,", NAN_REPLACEMENT);
            }
            else {
                if (CMath::is_finite(get_p(i)))
                    fprintf(file, "%e", (double)get_p(i));
                else
                    fprintf(file, "%f", NAN_REPLACEMENT);
            }
        }

        fprintf(file,"];\n\nq=[");
        for (i=0; i<N; i++)
        {
            if (i<N-1) {
                if (CMath::is_finite(get_q(i)))
                    fprintf(file, "%e,", (double)get_q(i));
                else
                    fprintf(file, "%f,", NAN_REPLACEMENT);
            }
            else {
                if (CMath::is_finite(get_q(i)))
                    fprintf(file, "%e", (double)get_q(i));
                else
                    fprintf(file, "%f", NAN_REPLACEMENT);
            }
        }
        fprintf(file,"];\n\na=[");

        for (i=0; i<N; i++)
        {
            fprintf(file, "\t[");

            for (j=0; j<N; j++)
            {
                if (j<N-1) {
                    if (CMath::is_finite(get_a(i,j)))
                        fprintf(file, "%e,", (double)get_a(i,j));
                    else
                        fprintf(file, "%f,", NAN_REPLACEMENT);
                }
                else {
                    if (CMath::is_finite(get_a(i,j)))
                        fprintf(file, "%e];\n", (double)get_a(i,j));
                    else
                        fprintf(file, "%f];\n", NAN_REPLACEMENT);
                }
            }
        }

        fprintf(file,"  ];\n\nb=[");

        for (i=0; i<N; i++)
        {
            fprintf(file, "\t[");

            for (j=0; j<M; j++)
            {
                if (j<M-1) {
                    if (CMath::is_finite(get_b(i,j)))
                        fprintf(file, "%e,",  (double)get_b(i,j));
                    else
                        fprintf(file, "%f,", NAN_REPLACEMENT);
                }
                else {
                    if (CMath::is_finite(get_b(i,j)))
                        fprintf(file, "%e];\n", (double)get_b(i,j));
                    else
                        fprintf(file, "%f];\n", NAN_REPLACEMENT);
                }
            }
        }
        result= (fprintf(file,"  ];\n") == 5);
    }

    return result;
}

T_STATES* CHMM::get_path(int32_t dim, float64_t& prob)
{
    T_STATES* result = NULL;

    prob = best_path(dim);
    result = SG_MALLOC(T_STATES, p_observations->get_vector_length(dim));

    for (int32_t i=0; i<p_observations->get_vector_length(dim); i++)
        result[i]=PATH(dim)[i];

    return result;
}

bool CHMM::save_path(FILE* file)
{
    bool result=false;

    if (file)
    {
      for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
        {
          if (dim%100==0)
        SG_PRINT("%i..", dim)
          float64_t prob = best_path(dim);
          fprintf(file,"%i. path probability:%e\nstate sequence:\n", dim, prob);
          for (int32_t i=0; i<p_observations->get_vector_length(dim)-1; i++)
        fprintf(file,"%d ", PATH(dim)[i]);
          fprintf(file,"%d", PATH(dim)[p_observations->get_vector_length(dim)-1]);
          fprintf(file,"\n\n") ;
        }
      SG_DONE()
      result=true;
    }

    return result;
}

bool CHMM::save_likelihood_bin(FILE* file)
{
    bool result=false;

    if (file)
    {
        for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
        {
            float32_t prob= (float32_t) model_probability(dim);
            fwrite(&prob, sizeof(float32_t), 1, file);
        }
        result=true;
    }

    return result;
}

bool CHMM::save_likelihood(FILE* file)
{
    bool result=false;

    if (file)
    {
        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");

        fprintf(file, "P=[");
        for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
            fprintf(file, "%e ", (double) model_probability(dim));

        fprintf(file,"];");
        result=true;
    }

    return result;
}

#define FLOATWRITE(file, value) { float32_t rrr=float32_t(value); fwrite(&rrr, sizeof(float32_t), 1, file); num_floats++;}

bool CHMM::save_model_bin(FILE* file)
{
    int32_t i,j,q, num_floats=0 ;
    if (!model)
    {
        if (file)
        {
            // write id
            FLOATWRITE(file, (float32_t)CMath::INFTY);
            FLOATWRITE(file, (float32_t) 1);

            //derivates log(dp),log(dq)
            for (i=0; i<N; i++)
                FLOATWRITE(file, get_p(i));
            SG_INFO("wrote %i parameters for p\n",N)

            for (i=0; i<N; i++)
                FLOATWRITE(file, get_q(i)) ;
            SG_INFO("wrote %i parameters for q\n",N)

            //derivates log(da),log(db)
            for (i=0; i<N; i++)
                for (j=0; j<N; j++)
                    FLOATWRITE(file, get_a(i,j));
            SG_INFO("wrote %i parameters for a\n",N*N)

            for (i=0; i<N; i++)
                for (j=0; j<M; j++)
                    FLOATWRITE(file, get_b(i,j));
            SG_INFO("wrote %i parameters for b\n",N*M)

            // write id
            FLOATWRITE(file, (float32_t)CMath::INFTY);
            FLOATWRITE(file, (float32_t) 3);

            // write number of parameters
            FLOATWRITE(file, (float32_t) N);
            FLOATWRITE(file, (float32_t) N);
            FLOATWRITE(file, (float32_t) N*N);
            FLOATWRITE(file, (float32_t) N*M);
            FLOATWRITE(file, (float32_t) N);
            FLOATWRITE(file, (float32_t) M);
        } ;
    }
    else
    {
        if (file)
        {
            int32_t num_p, num_q, num_a, num_b ;
            // write id
            FLOATWRITE(file, (float32_t)CMath::INFTY);
            FLOATWRITE(file, (float32_t) 2);

            for (i=0; model->get_learn_p(i)>=0; i++)
                FLOATWRITE(file, get_p(model->get_learn_p(i)));
            num_p=i ;
            SG_INFO("wrote %i parameters for p\n",num_p)

            for (i=0; model->get_learn_q(i)>=0; i++)
                FLOATWRITE(file, get_q(model->get_learn_q(i)));
            num_q=i ;
            SG_INFO("wrote %i parameters for q\n",num_q)

            //derivates log(da),log(db)
            for (q=0; model->get_learn_a(q,1)>=0; q++)
            {
                i=model->get_learn_a(q,0) ;
                j=model->get_learn_a(q,1) ;
                FLOATWRITE(file, (float32_t)i);
                FLOATWRITE(file, (float32_t)j);
                FLOATWRITE(file, get_a(i,j));
            } ;
            num_a=q ;
            SG_INFO("wrote %i parameters for a\n",num_a)

            for (q=0; model->get_learn_b(q,0)>=0; q++)
            {
                i=model->get_learn_b(q,0) ;
                j=model->get_learn_b(q,1) ;
                FLOATWRITE(file, (float32_t)i);
                FLOATWRITE(file, (float32_t)j);
                FLOATWRITE(file, get_b(i,j));
            } ;
            num_b=q ;
            SG_INFO("wrote %i parameters for b\n",num_b)

            // write id
            FLOATWRITE(file, (float32_t)CMath::INFTY);
            FLOATWRITE(file, (float32_t) 3);

            // write number of parameters
            FLOATWRITE(file, (float32_t) num_p);
            FLOATWRITE(file, (float32_t) num_q);
            FLOATWRITE(file, (float32_t) num_a);
            FLOATWRITE(file, (float32_t) num_b);
            FLOATWRITE(file, (float32_t) N);
            FLOATWRITE(file, (float32_t) M);
        } ;
    } ;
    return true ;
}

bool CHMM::save_path_derivatives(FILE* logfile)
{
    int32_t dim,i,j;

    if (logfile)
    {
        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");
        fprintf(logfile,"%% dPr[...]=[ [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
        fprintf(logfile,"%%            [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
        fprintf(logfile,"%%                            .............................                                \n");
        fprintf(logfile,"%%            [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_MM]\n");
        fprintf(logfile,"%%          ];\n%%\n\ndPr(log()) = [\n");
    }
    else
        return false ;

    for (dim=0; dim<p_observations->get_num_vectors(); dim++)
    {
        best_path(dim);

        fprintf(logfile, "[ ");

        //derivates dlogp,dlogq
        for (i=0; i<N; i++)
            fprintf(logfile,"%e, ", (double) path_derivative_p(i,dim) );

        for (i=0; i<N; i++)
            fprintf(logfile,"%e, ", (double) path_derivative_q(i,dim) );

        //derivates dloga,dlogb
        for (i=0; i<N; i++)
            for (j=0; j<N; j++)
                fprintf(logfile, "%e,", (double) path_derivative_a(i,j,dim) );

        for (i=0; i<N; i++)
            for (j=0; j<M; j++)
                fprintf(logfile, "%e,", (double) path_derivative_b(i,j,dim) );

        fseek(logfile,ftell(logfile)-1,SEEK_SET);
        fprintf(logfile, " ];\n");
    }

    fprintf(logfile, "];");

    return true ;
}

bool CHMM::save_path_derivatives_bin(FILE* logfile)
{
    bool result=false;
    int32_t dim,i,j,q;
    float64_t prob=0 ;
    int32_t num_floats=0 ;

    float64_t sum_prob=0.0 ;
    if (!model)
        SG_WARNING("No definitions loaded -- writing derivatives of all weights\n")
    else
        SG_INFO("writing derivatives of changed weights only\n")

    for (dim=0; dim<p_observations->get_num_vectors(); dim++)
    {
        if (dim%100==0)
        {
            SG_PRINT(".")

        } ;

        prob=best_path(dim);
        sum_prob+=prob ;

        if (!model)
        {
            if (logfile)
            {
                // write prob
                FLOATWRITE(logfile, prob);

                for (i=0; i<N; i++)
                    FLOATWRITE(logfile, path_derivative_p(i,dim));

                for (i=0; i<N; i++)
                    FLOATWRITE(logfile, path_derivative_q(i,dim));

                for (i=0; i<N; i++)
                    for (j=0; j<N; j++)
                        FLOATWRITE(logfile, path_derivative_a(i,j,dim));

                for (i=0; i<N; i++)
                    for (j=0; j<M; j++)
                        FLOATWRITE(logfile, path_derivative_b(i,j,dim));

            }
        }
        else
        {
            if (logfile)
            {
                // write prob
                FLOATWRITE(logfile, prob);

                for (i=0; model->get_learn_p(i)>=0; i++)
                    FLOATWRITE(logfile, path_derivative_p(model->get_learn_p(i),dim));

                for (i=0; model->get_learn_q(i)>=0; i++)
                    FLOATWRITE(logfile, path_derivative_q(model->get_learn_q(i),dim));

                for (q=0; model->get_learn_a(q,0)>=0; q++)
                {
                    i=model->get_learn_a(q,0) ;
                    j=model->get_learn_a(q,1) ;
                    FLOATWRITE(logfile, path_derivative_a(i,j,dim));
                } ;

                for (q=0; model->get_learn_b(q,0)>=0; q++)
                {
                    i=model->get_learn_b(q,0) ;
                    j=model->get_learn_b(q,1) ;
                    FLOATWRITE(logfile, path_derivative_b(i,j,dim));
                } ;
            }
        } ;
    }
    save_model_bin(logfile) ;

    result=true;
    SG_PRINT("\n")
    return result;
}

bool CHMM::save_model_derivatives_bin(FILE* file)
{
    bool result=false;
    int32_t dim,i,j,q ;
    int32_t num_floats=0 ;

    if (!model)
        SG_WARNING("No definitions loaded -- writing derivatives of all weights\n")
    else
        SG_INFO("writing derivatives of changed weights only\n")

#ifdef USE_HMMPARALLEL
    int32_t num_threads = parallel->get_num_threads();
    pthread_t *threads=SG_MALLOC(pthread_t, num_threads);
    S_DIM_THREAD_PARAM *params=SG_MALLOC(S_DIM_THREAD_PARAM, num_threads);

    if (p_observations->get_num_vectors()<num_threads)
        num_threads=p_observations->get_num_vectors();
#endif

    for (dim=0; dim<p_observations->get_num_vectors(); dim++)
    {
        if (dim%20==0)
        {
            SG_PRINT(".")

        } ;

#ifdef USE_HMMPARALLEL
        if (dim%num_threads==0)
        {
            for (i=0; i<num_threads; i++)
            {
                if (dim+i<p_observations->get_num_vectors())
                {
                    params[i].hmm=this ;
                    params[i].dim=dim+i ;
                    pthread_create(&threads[i], NULL, bw_dim_prefetch, (void*)&params[i]) ;
                }
            }

            for (i=0; i<num_threads; i++)
            {
                if (dim+i<p_observations->get_num_vectors())
                    pthread_join(threads[i], NULL);
            }
        }
#endif

        float64_t prob=model_probability(dim) ;
        if (!model)
        {
            if (file)
            {
                // write prob
                FLOATWRITE(file, prob);

                //derivates log(dp),log(dq)
                for (i=0; i<N; i++)
                    FLOATWRITE(file, model_derivative_p(i,dim));

                for (i=0; i<N; i++)
                    FLOATWRITE(file, model_derivative_q(i,dim));

                //derivates log(da),log(db)
                for (i=0; i<N; i++)
                    for (j=0; j<N; j++)
                        FLOATWRITE(file, model_derivative_a(i,j,dim));

                for (i=0; i<N; i++)
                    for (j=0; j<M; j++)
                        FLOATWRITE(file, model_derivative_b(i,j,dim));

                if (dim==0)
                    SG_INFO("Number of parameters (including posterior prob.): %i\n", num_floats)
            } ;
        }
        else
        {
            if (file)
            {
                // write prob
                FLOATWRITE(file, prob);

                for (i=0; model->get_learn_p(i)>=0; i++)
                    FLOATWRITE(file, model_derivative_p(model->get_learn_p(i),dim));

                for (i=0; model->get_learn_q(i)>=0; i++)
                    FLOATWRITE(file, model_derivative_q(model->get_learn_q(i),dim));

                //derivates log(da),log(db)
                for (q=0; model->get_learn_a(q,1)>=0; q++)
                {
                    i=model->get_learn_a(q,0) ;
                    j=model->get_learn_a(q,1) ;
                    FLOATWRITE(file, model_derivative_a(i,j,dim));
                } ;

                for (q=0; model->get_learn_b(q,0)>=0; q++)
                {
                    i=model->get_learn_b(q,0) ;
                    j=model->get_learn_b(q,1) ;
                    FLOATWRITE(file, model_derivative_b(i,j,dim));
                } ;
                if (dim==0)
                    SG_INFO("Number of parameters (including posterior prob.): %i\n", num_floats)
            } ;
        } ;
    }
    save_model_bin(file) ;

#ifdef USE_HMMPARALLEL
    SG_FREE(threads);
    SG_FREE(params);
#endif

    result=true;
    SG_PRINT("\n")
    return result;
}

bool CHMM::save_model_derivatives(FILE* file)
{
    bool result=false;
    int32_t dim,i,j;

    if (file)
    {

        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");
        fprintf(file,"%% dPr[...]=[ [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
        fprintf(file,"%%            [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_NN]\n");
        fprintf(file,"%%                            .............................                                \n");
        fprintf(file,"%%            [dp_1,...,dp_N,dq_1,...,dq_N, da_11,da_12,..,da_1N,..,da_NN, db_11,.., db_MM]\n");
        fprintf(file,"%%          ];\n%%\n\nlog(dPr) = [\n");


        for (dim=0; dim<p_observations->get_num_vectors(); dim++)
        {
            fprintf(file, "[ ");

            //derivates log(dp),log(dq)
            for (i=0; i<N; i++)
                fprintf(file,"%e, ", (double) model_derivative_p(i, dim) );     //log (dp)

            for (i=0; i<N; i++)
                fprintf(file,"%e, ", (double) model_derivative_q(i, dim) ); //log (dq)

            //derivates log(da),log(db)
            for (i=0; i<N; i++)
                for (j=0; j<N; j++)
                    fprintf(file, "%e,", (double) model_derivative_a(i,j,dim) );

            for (i=0; i<N; i++)
                for (j=0; j<M; j++)
                    fprintf(file, "%e,", (double) model_derivative_b(i,j,dim) );

            fseek(file,ftell(file)-1,SEEK_SET);
            fprintf(file, " ];\n");
        }


        fprintf(file, "];");

        result=true;
    }
    return result;
}

bool CHMM::check_model_derivatives_combined()
{
    //  bool result=false;
    const float64_t delta=5e-4 ;

    int32_t i ;
    //derivates log(da)
    /*  for (i=0; i<N; i++)
        {
        for (int32_t j=0; j<N; j++)
        {
        float64_t old_a=get_a(i,j) ;

        set_a(i,j, log(exp(old_a)-delta)) ;
        invalidate_model() ;
        float64_t prob_old=exp(model_probability(-1)*p_observations->get_num_vectors()) ;

        set_a(i,j, log(exp(old_a)+delta)) ;
        invalidate_model() ;
        float64_t prob_new=exp(model_probability(-1)*p_observations->get_num_vectors());

        float64_t deriv = (prob_new-prob_old)/(2*delta) ;

        set_a(i,j, old_a) ;
        invalidate_model() ;

        float64_t prod_prob=model_probability(-1)*p_observations->get_num_vectors() ;

        float64_t deriv_calc=0 ;
        for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
        deriv_calc+=exp(model_derivative_a(i, j, dim)+
        prod_prob-model_probability(dim)) ;

        SG_DEBUG("da(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc,  deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
        } ;
        } ;*/
    //derivates log(db)
    i=0;//for (i=0; i<N; i++)
    {
        for (int32_t j=0; j<M; j++)
        {
            float64_t old_b=get_b(i,j) ;

            set_b(i,j, log(exp(old_b)-delta)) ;
            invalidate_model() ;
            float64_t prob_old=(model_probability(-1)*p_observations->get_num_vectors()) ;

            set_b(i,j, log(exp(old_b)+delta)) ;
            invalidate_model() ;
            float64_t prob_new=(model_probability(-1)*p_observations->get_num_vectors());

            float64_t deriv = (prob_new-prob_old)/(2*delta) ;

            set_b(i,j, old_b) ;
            invalidate_model() ;

            float64_t deriv_calc=0 ;
            for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
            {
                deriv_calc+=exp(model_derivative_b(i, j, dim)-model_probability(dim)) ;
                if (j==1)
                    SG_INFO("deriv_calc[%i]=%e\n",dim,deriv_calc)
            } ;

            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)
        } ;
    } ;
    return true ;
}

bool CHMM::check_model_derivatives()
{
    bool result=false;
    const float64_t delta=3e-4 ;

    for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
    {
        int32_t i ;
        //derivates log(dp),log(dq)
        for (i=0; i<N; i++)
        {
            for (int32_t j=0; j<N; j++)
            {
                float64_t old_a=get_a(i,j) ;

                set_a(i,j, log(exp(old_a)-delta)) ;
                invalidate_model() ;
                float64_t prob_old=exp(model_probability(dim)) ;

                set_a(i,j, log(exp(old_a)+delta)) ;
                invalidate_model() ;
                float64_t prob_new=exp(model_probability(dim));

                float64_t deriv = (prob_new-prob_old)/(2*delta) ;

                set_a(i,j, old_a) ;
                invalidate_model() ;
                float64_t deriv_calc=exp(model_derivative_a(i, j, dim)) ;

                SG_DEBUG("da(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc,  deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
                invalidate_model() ;
            } ;
        } ;
        for (i=0; i<N; i++)
        {
            for (int32_t j=0; j<M; j++)
            {
                float64_t old_b=get_b(i,j) ;

                set_b(i,j, log(exp(old_b)-delta)) ;
                invalidate_model() ;
                float64_t prob_old=exp(model_probability(dim)) ;

                set_b(i,j, log(exp(old_b)+delta)) ;
                invalidate_model() ;
                float64_t prob_new=exp(model_probability(dim));

                float64_t deriv = (prob_new-prob_old)/(2*delta) ;

                set_b(i,j, old_b) ;
                invalidate_model() ;
                float64_t deriv_calc=exp(model_derivative_b(i, j, dim));

                SG_DEBUG("db(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/(deriv_calc))
            } ;
        } ;

#ifdef TEST
        for (i=0; i<N; i++)
        {
            float64_t old_p=get_p(i) ;

            set_p(i, log(exp(old_p)-delta)) ;
            invalidate_model() ;
            float64_t prob_old=exp(model_probability(dim)) ;

            set_p(i, log(exp(old_p)+delta)) ;
            invalidate_model() ;
            float64_t prob_new=exp(model_probability(dim));
            float64_t deriv = (prob_new-prob_old)/(2*delta) ;

            set_p(i, old_p) ;
            invalidate_model() ;
            float64_t deriv_calc=exp(model_derivative_p(i, dim));

            //if (fabs(deriv_calc_old-deriv)>1e-4)
            SG_DEBUG("dp(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
        } ;
        for (i=0; i<N; i++)
        {
            float64_t old_q=get_q(i) ;

            set_q(i, log(exp(old_q)-delta)) ;
            invalidate_model() ;
            float64_t prob_old=exp(model_probability(dim)) ;

            set_q(i, log(exp(old_q)+delta)) ;
            invalidate_model() ;
            float64_t prob_new=exp(model_probability(dim));

            float64_t deriv = (prob_new-prob_old)/(2*delta) ;

            set_q(i, old_q) ;
            invalidate_model() ;
            float64_t deriv_calc=exp(model_derivative_q(i, dim));

            //if (fabs(deriv_calc_old-deriv)>1e-4)
            SG_DEBUG("dq(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
        } ;
#endif
    }
    return result;
}

#ifdef USE_HMMDEBUG
bool CHMM::check_path_derivatives()
{
    bool result=false;
    const float64_t delta=1e-4 ;

    for (int32_t dim=0; dim<p_observations->get_num_vectors(); dim++)
    {
        int32_t i ;
        //derivates log(dp),log(dq)
        for (i=0; i<N; i++)
        {
            for (int32_t j=0; j<N; j++)
            {
                float64_t old_a=get_a(i,j) ;

                set_a(i,j, log(exp(old_a)-delta)) ;
                invalidate_model() ;
                float64_t prob_old=best_path(dim) ;

                set_a(i,j, log(exp(old_a)+delta)) ;
                invalidate_model() ;
                float64_t prob_new=best_path(dim);

                float64_t deriv = (prob_new-prob_old)/(2*delta) ;

                set_a(i,j, old_a) ;
                invalidate_model() ;
                float64_t deriv_calc=path_derivative_a(i, j, dim) ;

                SG_DEBUG("da(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc,  deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
            } ;
        } ;
        for (i=0; i<N; i++)
        {
            for (int32_t j=0; j<M; j++)
            {
                float64_t old_b=get_b(i,j) ;

                set_b(i,j, log(exp(old_b)-delta)) ;
                invalidate_model() ;
                float64_t prob_old=best_path(dim) ;

                set_b(i,j, log(exp(old_b)+delta)) ;
                invalidate_model() ;
                float64_t prob_new=best_path(dim);

                float64_t deriv = (prob_new-prob_old)/(2*delta) ;

                set_b(i,j, old_b) ;
                invalidate_model() ;
                float64_t deriv_calc=path_derivative_b(i, j, dim);

                SG_DEBUG("db(%i,%i) = %e:%e\t (%1.5f%%)\n", i,j, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/(deriv_calc))
            } ;
        } ;

        for (i=0; i<N; i++)
        {
            float64_t old_p=get_p(i) ;

            set_p(i, log(exp(old_p)-delta)) ;
            invalidate_model() ;
            float64_t prob_old=best_path(dim) ;

            set_p(i, log(exp(old_p)+delta)) ;
            invalidate_model() ;
            float64_t prob_new=best_path(dim);
            float64_t deriv = (prob_new-prob_old)/(2*delta) ;

            set_p(i, old_p) ;
            invalidate_model() ;
            float64_t deriv_calc=path_derivative_p(i, dim);

            //if (fabs(deriv_calc_old-deriv)>1e-4)
            SG_DEBUG("dp(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
        } ;
        for (i=0; i<N; i++)
        {
            float64_t old_q=get_q(i) ;

            set_q(i, log(exp(old_q)-delta)) ;
            invalidate_model() ;
            float64_t prob_old=best_path(dim) ;

            set_q(i, log(exp(old_q)+delta)) ;
            invalidate_model() ;
            float64_t prob_new=best_path(dim);

            float64_t deriv = (prob_new-prob_old)/(2*delta) ;

            set_q(i, old_q) ;
            invalidate_model() ;
            float64_t deriv_calc=path_derivative_q(i, dim);

            //if (fabs(deriv_calc_old-deriv)>1e-4)
            SG_DEBUG("dq(%i) = %e:%e\t (%1.5f%%)\n", i, deriv_calc, deriv, 100.0*(deriv-deriv_calc)/deriv_calc)
        } ;
    }
    return result;
}
#endif // USE_HMMDEBUG

//normalize model (sum to one constraint)
void CHMM::normalize(bool keep_dead_states)
{
    int32_t i,j;
    const float64_t INF=-1e10;
    float64_t sum_p =INF;

    for (i=0; i<N; i++)
    {
        sum_p=CMath::logarithmic_sum(sum_p, get_p(i));

        float64_t sum_b =INF;
        float64_t sum_a =get_q(i);

        for (j=0; j<N; j++)
            sum_a=CMath::logarithmic_sum(sum_a, get_a(i,j));

        if (sum_a>CMath::ALMOST_NEG_INFTY/N || (!keep_dead_states) )
        {
            for (j=0; j<N; j++)
                set_a(i,j, get_a(i,j)-sum_a);
            set_q(i, get_q(i)-sum_a);
        }

        for (j=0; j<M; j++)
            sum_b=CMath::logarithmic_sum(sum_b, get_b(i,j));
        for (j=0; j<M; j++)
            set_b(i,j, get_b(i,j)-sum_b);
    }

    for (i=0; i<N; i++)
        set_p(i, get_p(i)-sum_p);

    invalidate_model();
}

bool CHMM::append_model(CHMM* app_model)
{
    bool result=false;
    const int32_t num_states=app_model->get_N();
    int32_t i,j;

    SG_DEBUG("cur N:%d M:%d\n", N, M)
    SG_DEBUG("old N:%d M:%d\n", app_model->get_N(), app_model->get_M())
    if (app_model->get_M() == get_M())
    {
        float64_t* n_p=SG_MALLOC(float64_t, N+num_states);
        float64_t* n_q=SG_MALLOC(float64_t, N+num_states);
        float64_t* n_a=SG_MALLOC(float64_t, (N+num_states)*(N+num_states));
        //SG_PRINT("size n_b: %d\n", (N+num_states)*M)
        float64_t* n_b=SG_MALLOC(float64_t, (N+num_states)*M);

        //clear n_x
        for (i=0; i<N+num_states; i++)
        {
            n_p[i]=-CMath::INFTY;
            n_q[i]=-CMath::INFTY;

            for (j=0; j<N+num_states; j++)
                n_a[(N+num_states)*i+j]=-CMath::INFTY;

            for (j=0; j<M; j++)
                n_b[M*i+j]=-CMath::INFTY;
        }

        //copy models first
        // warning pay attention to the ordering of
        // transition_matrix_a, observation_matrix_b !!!

        // cur_model
        for (i=0; i<N; i++)
        {
            n_p[i]=get_p(i);

            for (j=0; j<N; j++)
                n_a[(N+num_states)*j+i]=get_a(i,j);

            for (j=0; j<M; j++)
            {
                n_b[M*i+j]=get_b(i,j);
            }
        }

        // append_model
        for (i=0; i<app_model->get_N(); i++)
        {
            n_q[i+N]=app_model->get_q(i);

            for (j=0; j<app_model->get_N(); j++)
                n_a[(N+num_states)*(j+N)+(i+N)]=app_model->get_a(i,j);
            for (j=0; j<app_model->get_M(); j++)
                n_b[M*(i+N)+j]=app_model->get_b(i,j);
        }


        // transition to the two and back
        for (i=0; i<N; i++)
        {
            for (j=N; j<N+num_states; j++)
                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]);
        }

        free_state_dependend_arrays();
        N+=num_states;

        alloc_state_dependend_arrays();

        //delete + adjust pointers
        SG_FREE(initial_state_distribution_p);
        SG_FREE(end_state_distribution_q);
        SG_FREE(transition_matrix_a);
        SG_FREE(observation_matrix_b);

        transition_matrix_a=n_a;
        observation_matrix_b=n_b;
        initial_state_distribution_p=n_p;
        end_state_distribution_q=n_q;

        SG_WARNING("not normalizing anymore, call normalize_hmm to make sure the hmm is valid!!\n")
        invalidate_model();
    }
    else
        SG_ERROR("number of observations is different for append model, doing nothing!\n")

    return result;
}

bool CHMM::append_model(CHMM* app_model, float64_t* cur_out, float64_t* app_out)
{
    bool result=false;
    const int32_t num_states=app_model->get_N()+2;
    int32_t i,j;

    if (app_model->get_M() == get_M())
    {
        float64_t* n_p=SG_MALLOC(float64_t, N+num_states);
        float64_t* n_q=SG_MALLOC(float64_t, N+num_states);
        float64_t* n_a=SG_MALLOC(float64_t, (N+num_states)*(N+num_states));
        //SG_PRINT("size n_b: %d\n", (N+num_states)*M)
        float64_t* n_b=SG_MALLOC(float64_t, (N+num_states)*M);

        //clear n_x
        for (i=0; i<N+num_states; i++)
        {
            n_p[i]=-CMath::INFTY;
            n_q[i]=-CMath::INFTY;

            for (j=0; j<N+num_states; j++)
                n_a[(N+num_states)*j+i]=-CMath::INFTY;

            for (j=0; j<M; j++)
                n_b[M*i+j]=-CMath::INFTY;
        }

        //copy models first
        // warning pay attention to the ordering of
        // transition_matrix_a, observation_matrix_b !!!

        // cur_model
        for (i=0; i<N; i++)
        {
            n_p[i]=get_p(i);

            for (j=0; j<N; j++)
                n_a[(N+num_states)*j+i]=get_a(i,j);

            for (j=0; j<M; j++)
            {
                n_b[M*i+j]=get_b(i,j);
            }
        }

        // append_model
        for (i=0; i<app_model->get_N(); i++)
        {
            n_q[i+N+2]=app_model->get_q(i);

            for (j=0; j<app_model->get_N(); j++)
                n_a[(N+num_states)*(j+N+2)+(i+N+2)]=app_model->get_a(i,j);
            for (j=0; j<app_model->get_M(); j++)
                n_b[M*(i+N+2)+j]=app_model->get_b(i,j);
        }

        //initialize the two special states

        // output
        for (i=0; i<M; i++)
        {
            n_b[M*N+i]=cur_out[i];
            n_b[M*(N+1)+i]=app_out[i];
        }

        // transition to the two and back
        for (i=0; i<N+num_states; i++)
        {
            // the first state is only connected to the second
            if (i==N+1)
                n_a[(N+num_states)*i + N]=0;

            // only states of the cur_model can reach the
            // first state
            if ((uint32_t)i<(uint32_t)N)
                n_a[(N+num_states)*N+i]=get_q(i);

            // the second state is only connected to states of
            // the append_model (with probab app->p(i))
            if (i>=N+2)
                n_a[(N+num_states)*i+(N+1)]=app_model->get_p(i-(N+2));
        }

        free_state_dependend_arrays();
        N+=num_states;

        alloc_state_dependend_arrays();

        //delete + adjust pointers
        SG_FREE(initial_state_distribution_p);
        SG_FREE(end_state_distribution_q);
        SG_FREE(transition_matrix_a);
        SG_FREE(observation_matrix_b);

        transition_matrix_a=n_a;
        observation_matrix_b=n_b;
        initial_state_distribution_p=n_p;
        end_state_distribution_q=n_q;

        SG_WARNING("not normalizing anymore, call normalize_hmm to make sure the hmm is valid!!\n")
        invalidate_model();
    }

    return result;
}


void CHMM::add_states(int32_t num_states, float64_t default_value)
{
    int32_t i,j;
    const float64_t MIN_RAND=1e-2; //this is the range of the random values for the new variables
    const float64_t MAX_RAND=2e-1;

    float64_t* n_p=SG_MALLOC(float64_t, N+num_states);
    float64_t* n_q=SG_MALLOC(float64_t, N+num_states);
    float64_t* n_a=SG_MALLOC(float64_t, (N+num_states)*(N+num_states));
    //SG_PRINT("size n_b: %d\n", (N+num_states)*M)
    float64_t* n_b=SG_MALLOC(float64_t, (N+num_states)*M);

    // warning pay attention to the ordering of
    // transition_matrix_a, observation_matrix_b !!!
    for (i=0; i<N; i++)
    {
        n_p[i]=get_p(i);
        n_q[i]=get_q(i);

        for (j=0; j<N; j++)
            n_a[(N+num_states)*j+i]=get_a(i,j);

        for (j=0; j<M; j++)
            n_b[M*i+j]=get_b(i,j);
    }

    for (i=N; i<N+num_states; i++)
    {
        n_p[i]=VAL_MACRO;
        n_q[i]=VAL_MACRO;

        for (j=0; j<N; j++)
            n_a[(N+num_states)*i+j]=VAL_MACRO;

        for (j=0; j<N+num_states; j++)
            n_a[(N+num_states)*j+i]=VAL_MACRO;

        for (j=0; j<M; j++)
            n_b[M*i+j]=VAL_MACRO;
    }
    free_state_dependend_arrays();
    N+=num_states;

    alloc_state_dependend_arrays();

    //delete + adjust pointers
    SG_FREE(initial_state_distribution_p);
    SG_FREE(end_state_distribution_q);
    SG_FREE(transition_matrix_a);
    SG_FREE(observation_matrix_b);

    transition_matrix_a=n_a;
    observation_matrix_b=n_b;
    initial_state_distribution_p=n_p;
    end_state_distribution_q=n_q;

    invalidate_model();
    normalize();
}

void CHMM::chop(float64_t value)
{
    for (int32_t i=0; i<N; i++)
    {
        int32_t j;

        if (exp(get_p(i)) < value)
            set_p(i, CMath::ALMOST_NEG_INFTY);

        if (exp(get_q(i)) < value)
            set_q(i, CMath::ALMOST_NEG_INFTY);

        for (j=0; j<N; j++)
        {
            if (exp(get_a(i,j)) < value)
                set_a(i,j, CMath::ALMOST_NEG_INFTY);
        }

        for (j=0; j<M; j++)
        {
            if (exp(get_b(i,j)) < value)
                set_b(i,j, CMath::ALMOST_NEG_INFTY);
        }
    }
    normalize();
    invalidate_model();
}

bool CHMM::linear_train(bool right_align)
{
    if (p_observations)
    {
        int32_t histsize=(get_M()*get_N());
        int32_t* hist=SG_MALLOC(int32_t, histsize);
        int32_t* startendhist=SG_MALLOC(int32_t, get_N());
        int32_t i,dim;

        ASSERT(p_observations->get_max_vector_length()<=get_N())

        for (i=0; i<histsize; i++)
            hist[i]=0;

        for (i=0; i<get_N(); i++)
            startendhist[i]=0;

        if (right_align)
        {
            for (dim=0; dim<p_observations->get_num_vectors(); dim++)
            {
                int32_t len=0;
                bool free_vec;
                uint16_t* obs=p_observations->get_feature_vector(dim, len, free_vec);

                ASSERT(len<=get_N())
                startendhist[(get_N()-len)]++;

                for (i=0;i<len;i++)
                    hist[(get_N()-len+i)*get_M() + *obs++]++;

                p_observations->free_feature_vector(obs, dim, free_vec);
            }

            set_q(get_N()-1, 1);
            for (i=0; i<get_N()-1; i++)
                set_q(i, 0);

            for (i=0; i<get_N(); i++)
                set_p(i, startendhist[i]+PSEUDO);

            for (i=0;i<get_N();i++)
            {
                for (int32_t j=0; j<get_N(); j++)
                {
                    if (i==j-1)
                        set_a(i,j, 1);
                    else
                        set_a(i,j, 0);
                }
            }
        }
        else
        {
            for (dim=0; dim<p_observations->get_num_vectors(); dim++)
            {
                int32_t len=0;
                bool free_vec;
                uint16_t* obs=p_observations->get_feature_vector(dim, len, free_vec);

                ASSERT(len<=get_N())
                for (i=0;i<len;i++)
                    hist[i*get_M() + *obs++]++;

                startendhist[len-1]++;

                p_observations->free_feature_vector(obs, dim, free_vec);
            }

            set_p(0, 1);
            for (i=1; i<get_N(); i++)
                set_p(i, 0);

            for (i=0; i<get_N(); i++)
                set_q(i, startendhist[i]+PSEUDO);

            int32_t total=p_observations->get_num_vectors();

            for (i=0;i<get_N();i++)
            {
                total-= startendhist[i] ;

                for (int32_t j=0; j<get_N(); j++)
                {
                    if (i==j-1)
                        set_a(i,j, total+PSEUDO);
                    else
                        set_a(i,j, 0);
                }
            }
            ASSERT(total==0)
        }

        for (i=0;i<get_N();i++)
        {
            for (int32_t j=0; j<get_M(); j++)
            {
                float64_t sum=0;
                int32_t offs=i*get_M()+ p_observations->get_masked_symbols((uint16_t) j, (uint8_t) 254);

                for (int32_t k=0; k<p_observations->get_original_num_symbols(); k++)
                    sum+=hist[offs+k];

                set_b(i,j, (PSEUDO+hist[i*get_M()+j])/(sum+PSEUDO*p_observations->get_original_num_symbols()));
            }
        }

        SG_FREE(hist);
        SG_FREE(startendhist);
        convert_to_log();
        invalidate_model();
        return true;
    }
    else
        return false;
}

void CHMM::set_observation_nocache(CStringFeatures<uint16_t>* obs)
{
    ASSERT(obs)
    p_observations=obs;
    SG_REF(obs);

    if (obs)
        if (obs->get_num_symbols() > M)
            SG_ERROR("number of symbols in observation (%ld) larger than M (%d)\n", (long) obs->get_num_symbols(), M)

    if (!reused_caches)
    {
#ifdef USE_HMMPARALLEL_STRUCTURES
        for (int32_t i=0; i<parallel->get_num_threads(); i++)
        {
            SG_FREE(alpha_cache[i].table);
            SG_FREE(beta_cache[i].table);
            SG_FREE(states_per_observation_psi[i]);
            SG_FREE(path[i]);

            alpha_cache[i].table=NULL;
            beta_cache[i].table=NULL;
            states_per_observation_psi[i]=NULL;
            path[i]=NULL;
        } ;
#else
        SG_FREE(alpha_cache.table);
        SG_FREE(beta_cache.table);
        SG_FREE(states_per_observation_psi);
        SG_FREE(path);

        alpha_cache.table=NULL;
        beta_cache.table=NULL;
        states_per_observation_psi=NULL;
        path=NULL;

#endif //USE_HMMPARALLEL_STRUCTURES
    }

    invalidate_model();
}

void CHMM::set_observations(CStringFeatures<uint16_t>* obs, CHMM* lambda)
{
    ASSERT(obs)
    SG_REF(obs);
    p_observations=obs;

    /* from Distribution, necessary for calls to base class methods, like
     * get_log_likelihood_sample():
     */
    SG_REF(obs);
    features=obs;

    SG_DEBUG("num symbols alphabet: %ld\n", obs->get_alphabet()->get_num_symbols())
    SG_DEBUG("num symbols: %ld\n", obs->get_num_symbols())
    SG_DEBUG("M: %d\n", M)

    if (obs)
    {
        if (obs->get_num_symbols() > M)
        {
            SG_ERROR("number of symbols in observation (%ld) larger than M (%d)\n", (long) obs->get_num_symbols(), M)
        }
    }

    if (!reused_caches)
    {
#ifdef USE_HMMPARALLEL_STRUCTURES
        for (int32_t i=0; i<parallel->get_num_threads(); i++)
        {
            SG_FREE(alpha_cache[i].table);
            SG_FREE(beta_cache[i].table);
            SG_FREE(states_per_observation_psi[i]);
            SG_FREE(path[i]);

            alpha_cache[i].table=NULL;
            beta_cache[i].table=NULL;
            states_per_observation_psi[i]=NULL;
            path[i]=NULL;
        } ;
#else
        SG_FREE(alpha_cache.table);
        SG_FREE(beta_cache.table);
        SG_FREE(states_per_observation_psi);
        SG_FREE(path);

        alpha_cache.table=NULL;
        beta_cache.table=NULL;
        states_per_observation_psi=NULL;
        path=NULL;

#endif //USE_HMMPARALLEL_STRUCTURES
    }

    if (obs!=NULL)
    {
        int32_t max_T=obs->get_max_vector_length();

        if (lambda)
        {
#ifdef USE_HMMPARALLEL_STRUCTURES
            for (int32_t i=0; i<parallel->get_num_threads(); i++)
            {
                this->alpha_cache[i].table= lambda->alpha_cache[i].table;
                this->beta_cache[i].table=  lambda->beta_cache[i].table;
                this->states_per_observation_psi[i]=lambda->states_per_observation_psi[i] ;
                this->path[i]=lambda->path[i];
            } ;
#else
            this->alpha_cache.table= lambda->alpha_cache.table;
            this->beta_cache.table= lambda->beta_cache.table;
            this->states_per_observation_psi= lambda->states_per_observation_psi;
            this->path=lambda->path;
#endif //USE_HMMPARALLEL_STRUCTURES

            this->reused_caches=true;
        }
        else
        {
            this->reused_caches=false;
#ifdef USE_HMMPARALLEL_STRUCTURES
            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)
            for (int32_t i=0; i<parallel->get_num_threads(); i++)
            {
                if ((states_per_observation_psi[i]=SG_MALLOC(T_STATES,max_T*N))!=NULL)
                    SG_DEBUG("path_table[%i] successfully allocated\n",i)
                else
                    SG_ERROR("failed allocating memory for path_table[%i].\n",i)
                path[i]=SG_MALLOC(T_STATES, max_T);
            }
#else // no USE_HMMPARALLEL_STRUCTURES
            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)
            if ((states_per_observation_psi=SG_MALLOC(T_STATES,max_T*N)) != NULL)
                SG_DONE()
            else
                SG_ERROR("failed.\n")

            path=SG_MALLOC(T_STATES, max_T);
#endif // USE_HMMPARALLEL_STRUCTURES
#ifdef USE_HMMCACHE
            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)

#ifdef USE_HMMPARALLEL_STRUCTURES
            for (int32_t i=0; i<parallel->get_num_threads(); i++)
            {
                if ((alpha_cache[i].table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N))!=NULL)
                    SG_DEBUG("alpha_cache[%i].table successfully allocated\n",i)
                else
                    SG_ERROR("allocation of alpha_cache[%i].table failed\n",i)

                if ((beta_cache[i].table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N)) != NULL)
                    SG_DEBUG("beta_cache[%i].table successfully allocated\n",i)
                else
                    SG_ERROR("allocation of beta_cache[%i].table failed\n",i)
            } ;
#else // USE_HMMPARALLEL_STRUCTURES
            if ((alpha_cache.table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N)) != NULL)
                SG_DEBUG("alpha_cache.table successfully allocated\n")
            else
                SG_ERROR("allocation of alpha_cache.table failed\n")

            if ((beta_cache.table=SG_MALLOC(T_ALPHA_BETA_TABLE, max_T*N)) != NULL)
                SG_DEBUG("beta_cache.table successfully allocated\n")
            else
                SG_ERROR("allocation of beta_cache.table failed\n")

#endif // USE_HMMPARALLEL_STRUCTURES
#else // USE_HMMCACHE
#ifdef USE_HMMPARALLEL_STRUCTURES
            for (int32_t i=0; i<parallel->get_num_threads(); i++)
            {
                alpha_cache[i].table=NULL ;
                beta_cache[i].table=NULL ;
            } ;
#else //USE_HMMPARALLEL_STRUCTURES
            alpha_cache.table=NULL ;
            beta_cache.table=NULL ;
#endif //USE_HMMPARALLEL_STRUCTURES
#endif //USE_HMMCACHE
        }
    }

    //initialize pat/mod_prob as not calculated
    invalidate_model();
}

bool CHMM::permutation_entropy(int32_t window_width, int32_t sequence_number)
{
    if (p_observations && window_width>0 &&
            ( sequence_number<0 || sequence_number < p_observations->get_num_vectors()))
    {
        int32_t min_sequence=sequence_number;
        int32_t max_sequence=sequence_number;

        if (sequence_number<0)
        {
            min_sequence=0;
            max_sequence=p_observations->get_num_vectors();
            SG_INFO("numseq: %d\n", max_sequence)
        }

        SG_INFO("min_sequence: %d max_sequence: %d\n", min_sequence, max_sequence)
        for (sequence_number=min_sequence; sequence_number<max_sequence; sequence_number++)
        {
            int32_t sequence_length=0;
            bool free_vec;
            uint16_t* obs=p_observations->get_feature_vector(sequence_number, sequence_length, free_vec);

            int32_t histsize=get_M();
            int64_t* hist=SG_MALLOC(int64_t, histsize);
            int32_t i,j;

            for (i=0; i<sequence_length-window_width; i++)
            {
                for (j=0; j<histsize; j++)
                    hist[j]=0;

                uint16_t* ptr=&obs[i];
                for (j=0; j<window_width; j++)
                {
                    hist[*ptr++]++;
                }

                float64_t perm_entropy=0;
                for (j=0; j<get_M(); j++)
                {
                    float64_t p=
                        (((float64_t) hist[j])+PSEUDO)/
                        (window_width+get_M()*PSEUDO);
                    perm_entropy+=p*log(p);
                }

                SG_PRINT("%f\n", perm_entropy)
            }
            p_observations->free_feature_vector(obs, sequence_number, free_vec);

            SG_FREE(hist);
        }
        return true;
    }
    else
        return false;
}

float64_t CHMM::get_log_derivative(int32_t num_param, int32_t num_example)
{
    if (num_param<N)
        return model_derivative_p(num_param, num_example);
    else if (num_param<2*N)
        return model_derivative_q(num_param-N, num_example);
    else if (num_param<N*(N+2))
    {
        int32_t k=num_param-2*N;
        int32_t i=k/N;
        int32_t j=k%N;
        return model_derivative_a(i,j, num_example);
    }
    else if (num_param<N*(N+2+M))
    {
        int32_t k=num_param-N*(N+2);
        int32_t i=k/M;
        int32_t j=k%M;
        return model_derivative_b(i,j, num_example);
    }

    ASSERT(false)
    return -1;
}

float64_t CHMM::get_log_model_parameter(int32_t num_param)
{
    if (num_param<N)
        return get_p(num_param);
    else if (num_param<2*N)
        return get_q(num_param-N);
    else if (num_param<N*(N+2))
        return transition_matrix_a[num_param-2*N];
    else if (num_param<N*(N+2+M))
        return observation_matrix_b[num_param-N*(N+2)];

    ASSERT(false)
    return -1;
}


//convergence criteria  -tobeadjusted-
bool CHMM::converged(float64_t x, float64_t y)
{
    float64_t diff=y-x;
    float64_t absdiff=fabs(diff);

    SG_INFO("\n #%03d\tbest result so far: %G (eps: %f)", iteration_count, y, diff)

    if (iteration_count--==0 || (absdiff<epsilon && conv_it<=0))
    {
        iteration_count=iterations;
        SG_INFO("...finished\n")
        conv_it=5;
        return true;
    }
    else
    {
        if (absdiff<epsilon)
            conv_it--;
        else
            conv_it=5;

        return false;
    }
}

bool CHMM::baum_welch_viterbi_train(BaumWelchViterbiType type)
{
    CHMM* estimate=new CHMM(this);
    CHMM* working=this;
    float64_t prob_max=-CMath::INFTY;
    float64_t prob=-CMath::INFTY;
    float64_t prob_train=CMath::ALMOST_NEG_INFTY;
    iteration_count=iterations;

    while (!converged(prob, prob_train) && (!CSignal::cancel_computations()))
    {
        CMath::swap(working, estimate);
        prob=prob_train;

        switch (type) {
            case BW_NORMAL:
                working->estimate_model_baum_welch(estimate); break;
            case BW_TRANS:
                working->estimate_model_baum_welch_trans(estimate); break;
            case BW_DEFINED:
                working->estimate_model_baum_welch_defined(estimate); break;
            case VIT_NORMAL:
                working->estimate_model_viterbi(estimate); break;
            case VIT_DEFINED:
                working->estimate_model_viterbi_defined(estimate); break;
        }
        prob_train=estimate->model_probability();

        if (prob_max<prob_train)
            prob_max=prob_train;
    }

    if (estimate == this)
    {
        estimate->copy_model(working);
        delete working;
    }
    else
        delete estimate;

    return true;
}