mleHMMintern <-
function(x,m,G,f,pn2pw,pw2pn,startParameters,startGamma,...){
    
  #global variables
  nState = length(pn2pw(startParameters[[1]]))
  
  #print(nState)
  
  #calculates stationary distribution
  stat_distr = function(mat){
    m0=dim(mat)[1]
    I_m = diag(rep(1,m0))
    U = matrix(rep(1,m0*m0),nrow=m0)
    M = I_m - mat + U
    One = matrix(rep(1,m0),nrow=1)
    erg = One %*% solve(M)
    return(as.vector(erg))
  }
  
  #given the restriction gamma = lambda_G(gamma)
  # the tpm gamma can be parametrized in nu and pList
  
  #transforms (k x k) tpm to (r x r) tpm and vector of probabilities
  gamma2nuP = function(gamma,G){
    
    #number of sets in partition
    r = length(G)  
    
    #matrix of group transitions
    nu = NA * matrix(0,ncol=r,nrow=r)
    #list of inner group probabilities
    pList = list()
    
    #calculates the probability from group l -> m
    probOfGroupTransition = function(l,m,G,gamma){
      #stationary distribution
      p0 = stat_distr(gamma)
      
      #probability for G_l
      sumGroup_l = sum(p0[G[[l]]])
      
      erg=0
      for(g in G[[m]]){
        for(h in G[[l]]){
          erg = erg + p0[h]/sumGroup_l * gamma[h,g]  
        } 
      }
      return(erg)
    }
    
    
    #set starting parameter for nu matrix
    for(i in 1:r){
      for(j in 1:r){
        nu[i,j] = probOfGroupTransition(l=i,m=j,G=G,gamma=gamma)
      }
    }
    
    #set starting parameter for pList
    p0 = stat_distr(gamma)
    for(l in 1:r){
      pList[[l]] = numeric(length(G[[l]]))
      sumGroup_l = sum(p0[G[[l]]])
      pList[[l]] = p0[G[[l]]]/sumGroup_l
    }
    
    return(list(nu=nu,pList=pList))
    
  }

  #transforms (r x r) tpm and vector of probabilities to (k x k) tpm 
  nuP2Gamma = function(nuPList,G){
    
    nu = nuPList$nu
    pList = nuPList$pList
    m=length(unlist(pList))
    
    getGroupOfIndex = function(i,G){
      r=length(G)
      for(l in 1:r){
        if(any(G[[l]]==i)) return(l)  
      }
    }
    
    getPositionInGroup = function(i,G){  
      group = getGroupOfIndex(i,G)
      return(which(G[[group]]==i))  
    }
    
    
    #set gamma matrix back 
    gamma = NA * matrix(0,ncol=m,nrow=m)
    
    for(i in 1:m){
      for(j in 1:m){
        gamma[i,j] = nu[getGroupOfIndex(i,G),getGroupOfIndex(j,G)] * pList[[getGroupOfIndex(j,G)]][getPositionInGroup(j,G)]  
      }
    }
    
    return(gamma)
    
  }
  
  #transformations to transform from/to restricted/unrestricted parameter space
  # the transformations for tpms and p's are allways the same
  # the transformations for the state dep. parameter depends on the family of the distribution
  
  ####
  #transforms tpm to vector
  tpmPn2pw = function(mat){
    m0  = dim(mat)[1]
    if(m0>1){
      foo = log(mat/diag(mat))
      tmat = as.vector(foo[!diag(m0)])
      return(tmat)  
    }
    else
      return(NULL)
  }
  
  #transforms vector to tpm
  tpmPw2pn = function(param,m0){
    epar = exp(param)
    mat = diag(m0)
    mat[!mat] = epar[1:(m0*m0-m0)]
    mat = mat/apply(mat,1,sum)
    mat
  }
  ####
  
  
  ####
  #transforms list of probability vectors to vector
  probPn2pw = function(pList,G){
    r = length(G)
    param = NULL
    
    for(l in 1:r){
      prob=pList[[l]]
      if(length(prob>1))
        param = c(param, log( prob[-1] / ( 1-sum(prob[-1]) ) ) )
      else
        param = c(param,NULL) 
    }
    return(param)
  }
  
  #transforms vector to list os probability vectors
  probPw2pn = function(param,G){
    pList = list()
    r = length(G)
    epar = exp(param)
    for(l in 1:r){
      if(length(G[[l]])==1)
        pList[[l]]=1
      else{
        nTMP = length(G[[l]])-1 #the first nTMP values of epar are relevant for the present group
        eParTMP = epar[1:nTMP]
        pTMP = eParTMP/(1+sum(eParTMP))
        pList[[l]]=c(1-sum(pTMP),pTMP) 
        epar = epar[-(1:nTMP)]
      }
    }
    return(pList)
  }
  ####
  
  ####
  #transforms all natural state dep. parameters to working parameters
  statePn2pw = function(parameters){
    stateWorking = NULL
    m=length(parameters)
    for(i in 1:m){
      stateWorking = c(stateWorking,pn2pw(parameters[[i]]))  
    }
    return(stateWorking)
  }
  
  #transforms all statedep. working paramaters to natural parameters
  statePw2pn = function(stateWorking){
    
    parameters = list()
    
    for(ind in 1:m){

      parameters[[ind]] = pw2pn(stateWorking[1:nState])
      stateWorking = stateWorking[-(1:nState)]
    }
    
    return(parameters)
  }
  ####
  
  ####
  #transformation of the complete parameter to a unconstranied vector
  allPn2Pw = function(nuPList,parameters,G){
  
    #read values
    r = dim(G)
    nu = nuPList$nu
    pList = nuPList$pList
    m = length(unlist(pList))
    
    
    #set working parameter
    #tpm nu
    workingTpm = tpmPn2pw(nu)
    #pList
    workingPlist = probPn2pw(pList=pList,G=G)
    #state dependent parameter
    workingState = statePn2pw(parameters=parameters)
    
    return(c(workingTpm,workingPlist,workingState))
    
  }
  

  allPw2pn = function(parameter,G){
    
    r = length(G)
    m = length(unlist(G))
    #tpm
    if(r==1){
      nu=matrix(1)
    }
    else{
      nu = tpmPw2pn(param=parameter[1:(r^2-r)],m0=r)
      parameter = parameter[-(1:(r^2-r))]
    }
    

    
    #pList
    z = 0; for(l in 1:r){ nG = length(G[[l]]);if(nG!=1) z=z+nG-1 }
    nParPlist = z

    if(z>0){  
      pList = probPw2pn(param=parameter[1:nParPlist],G=G)
      parameter = parameter[-(1:nParPlist)]
    }
    else{
      pList = probPw2pn(param=0,G=G)  
    }
    
    nuPList = list(nu=nu,pList=pList)
    
    stateParameters = statePw2pn(stateWorking=parameter)
    
    return(list(nuPList=nuPList,stateParameters=stateParameters))
    
  }
  ####
  

  
  norm.HMM.mllk <- function(parvect,x,m,G,...)
  {
    if(m==1) return(-Inf)
    
    #transform working parameters 2 natural parameters
    pn = allPw2pn(parameter=parvect,G=G)
    gamma = nuP2Gamma(nuPList = pn$nuPList, G = G)
    stateParameters = pn$stateParameters
    
    #stationary distribution
    try(delta <- stat_distr(gamma))
    if(class(try(delta <- stat_distr(gamma))) == "try-error")
    delta=rep(1/m,m)
    #print(delta)
    
    #no. of observations
    n = ifelse(test=is.matrix(x),yes=n<-dim(x)[1],no=n<-length(x))
    
    
    #state dependent density values
    allprobs= matrix(NA,ncol=m,nrow=n)
    for(j in 1:m){
      allprobs[,j] = f(x=x,par=stateParameters[[j]])
    }
    
    allprobs   <- ifelse(!is.na(allprobs),allprobs,1)
    lscale     <- 0
    foo        <- delta
    for (i in 1:n)
    {
      foo    <- foo%*%gamma*allprobs[i,]
      sumfoo <- sum(foo)
      lscale <- lscale+log(sumfoo)
      foo    <- foo/sumfoo
    }
    mllk       <- -(lscale)
    mllk
  }  
  
  #print("OK1")
  
  #print(statePn2pw(startParameters))
  
  parvect0 = allPn2Pw(nuPList=gamma2nuP(gamma=startGamma,G=G),parameters=startParameters,G=G)
  
  #print("OK2")
  
  mod = nlm(norm.HMM.mllk,parvect0,x=x,m=m,G=G,...)
  pn        <- allPw2pn(parameter=mod$estimate,G=G)
  mllk      <- mod$minimum
  np        <- length(parvect0)
  AIC       <- 2*(mllk+np)
  n         <- sum(!is.na(x))
  BIC       <- 2*mllk+np*log(n)
  par_est = list(gamma  = nuP2Gamma(nuPList=pn$nuPList,G=G), stateDepParameters = pn$stateParameters)
  
  
  return(list(mle=par_est,code=mod$code,mllk=mllk,AIC=AIC,BIC=BIC))
    
}
