viterbiHMM <-
function(x, hmm) {                                                           
  
  G = hmm@internVariables$G
  stateParameters = hmm@stateParameters

  nuPList = list(nu = hmm@transitionMatrix, pList = hmm@probList)
  
  gamma = nuP2Gamma(nuPList,G=hmm@internVariables$G)

  distribution = hmm@internVariables$distribution
  

  
  m=dim(gamma)[1]
  
  #reparametrisation
  stateParameters = extern2intern(stateParameters,distribution,m)
  
  #if(missing(G)) G = as.list(1:m)
  
  #no. of observations
  n = ifelse(test=is.matrix(x),yes=n<-dim(x)[1],no=n<-length(x))
  
  #old way -> no mixtures as state dep. distr. possible
  #allprobs <- matrix(NA,ncol=m,nrow=n)
  #for(i in 1:m) allprobs[,i] = f0(x,stateParameters[[i]])
  
  
  #mixtures for merged states
  dmixture = function(x,pList,G,stateParameters,l,distribution){
    
    states = G[[l]]
    probs = pList[[l]]
    
    np = length(probs)
    
    if(distribution=="poisson"){
      #density
      f0 = function(x,par){lambda=par;return(dpois(x,lambda=lambda))} 
    }
    
    if(distribution=="normal"){
      #density
      f0 = function(x,par){
        mu    = par[1]
        sigma = par[2]
        return(dnorm(x,mean=mu,sd=sigma))
      }
    }
    
    if(distribution=="mv_normal"){
      #density
      f0 = function(x,par){
        mu    = par$mu
        sigma = par$sigma
        return(dmnorm(x=x,mean=mu,varcov=sigma))
      }
    }
    
    erg=0
    for(i in 1:np){
      erg = erg + probs[i] * f0(x,stateParameters[[(states[i])]])
    }
    
    return(erg)
    
  }
  
  #calculates the state dep. mixture values given kxk tpm and partition G
  calculateStatedepDens  = function(x, stateParameters, gamma, G, distribution){
    
    nuP = gamma2nuP(gamma,G)
    nu = nuP$nu
    pList = nuP$pList
    
    #no. of states after merge, i.e. sets in G
    nNu = dim(nu)[1]
    
    
    #allprobs = matrix(NA,ncol=m,nrow=n)
    allprobs = dmixture(x,pList,G,stateParameters,l=1,distribution)
    
    if(m==1) return(allprobs)
    for(j in 2:nNu) allprobs = cbind(allprobs,dmixture(x,pList,G,stateParameters,l=j,distribution))
    
    return(allprobs)
  }
  
  #calculate state dep. density values
  allprobs = calculateStatedepDens(x, stateParameters, gamma, G, distribution)
  
  
  #tpm of merged model is nu
  nuP = gamma2nuP(gamma,G)
  nu = nuP$nu; m2 = dim(nu)[1]
  pList = nuP$pList
  
  delta<-solve(t(diag(m2)-nu+1),rep(1,m2))
  
  xi        <- matrix(0,n,m2)                                
  foo       <- delta*allprobs[1,]                          
  xi[1,]    <- foo/sum(foo)                                 
  for (i in 2:n)                                             
  {                                                        
    foo    <- apply(xi[i-1,]*nu,2,max)*allprobs[i,]      
    xi[i,] <- foo/sum(foo)                                   
  }                                                        
  iv<-numeric(n)                                             
  iv[n]     <-which.max(xi[n,])                             
  for (i in (n-1):1)                                         
    iv[i] <- which.max(nu[,iv[i+1]]*xi[i,])               
  iv                                                         
}
