findIndependencePartition <-
function(m,x,alpha=0.05,estimateIn,details = FALSE,hmm=NULL,...){

  
  if(class(x)=="numeric"){
    distribution="normal"
    mDim = 1
  }
  
  if(class(x)=="matrix"){
    if(dim(x)[2]==1){
      distribution="normal"   
      mDim = 1
    }
    else{
      distribution="mv_normal" 
      mDim = dim(x)[2]
    }
  }
  
 
  #initial partition
  Gstart  = as.list(1:m)
  
  #check if unrestricted estimateIn is provided
  if(missing(estimateIn)){
    cat("Estimate model with", m,"states.\n")
    mleFullModel = mleHMM(x=x,m=m,G=Gstart,hmm=hmm,...)
    hmm0 = mleFullModel$mle
    stateParam0 = hmm0@stateParameters
    gamma0 = hmm0@transitionMatrix
  }
  else{
    mleFullModel = estimateIn
    hmm0 = estimateIn$mle
    stateParam0 = hmm0@stateParameters
    gamma0 = hmm0@transitionMatrix
  }
  
  
    
  #given a partition G this functions finds all partitions combined to of the sets
  getAllSmallerPartitions = function(G){
    
    #set length of partition
    r= length(G)
    
    #get all possible combinations
    com = subsets(r, 2, v = 1:r)
    N=dim(com)[1]
    
    Glist = list()
    for(i in 1:N){
      j0 = com[i,] # j0 is vector, these groups should be combined
      Gtmp = G
      Gtmp[[(j0[1])]] = c(G[[(j0[1])]],G[[(j0[2])]])
      Gtmp[[(j0[2])]] = NULL
      Glist[[i]] = Gtmp
    }
    return(Glist)
  }
  

  #list of partitions during iteration
  Gpath = list(); Gpath[[m]] = Gstart
  
     
  #track path of p-Values
  pValuePath = NA * numeric(m)
  
  #track path of loglikelihod ratio's
  likelihoodRatioPath = NA *numeric(m)
  
  #which one is the final model?
  finalModel = m
  
  #track all models during reducing steps
  modelPath = list()
  modelPath[[m]] = mleFullModel
  
  
  #track path of log-likelihoods
  mllkPath = NA * numeric(m) 
  mllkPath[m] = mleFullModel$mllk
  
  estimateRestrictedModels = TRUE
  
  #for each partition the mle is calculated under the restriction
  #gamma = \lambda_G(gamma)
  if(estimateRestrictedModels == TRUE){
    
    #set a function which calculates all mle's of a partition list
    tmpMle = function(i,Glist,data,hmm,...){
      #print(gamma)
      Ghere = Glist[[i]]
      return(try(mleHMM(x=data,m=m,G=Glist[[i]],hmm=hmm,...)))
    } 
    
    for(l in (m-1):1){
      

      
      #last step: fit iid model
      if(l==1){
        cat("Estimate model with", l,"state.\n")
        testChar = class(try(modelPath[[l]] <- mleHMM(x=x,m=m,G=list(1:m),hmm=hmm0,...)))
        if(testChar=="try-error"){
          cat("Error in last estimation step (i.i.d. model).")
          break
        }
        Gpath[[l]] = list(1:m) 
        mllkPath[l] = modelPath[[l]]$mllk
        likelihoodRatioPath[l] = 2*(modelPath[[l]]$mllk - modelPath[[m]]$mllk)
        rNow = length(Gpath[[l]])     
        pValuePath[l] = pchisq(q=likelihoodRatioPath[l],
                               df=m^2-2*m-rNow^2+2*rNow,lower.tail=FALSE)
        #check if H_0: lambda_Gpath[[1]](gamma) = gamma can be rejected
        if(pValuePath[l]<alpha){
          finalModel = l+1
        }
        else
          finalModel = l  
        }
      #previous m-2 steps
      else{
        cat("Estimate models with", l,"states.\n")
        
        #calculate for all one-level-lower-partitions the mles under the restriction
        GlistTmp = getAllSmallerPartitions(G=Gpath[[l+1]])

        Btmp = length(GlistTmp)
        
        cl <- makeCluster(detectCores())

        #ANPASSEN AUF PACKGE, dann hier library(mergeHMM)
        clusterEvalQ(cl, {
          ## set up each worker.  Could also use clusterExport()
#           library(sn)
#           library(mvtnorm)
#           library(RHmm)
#           library(parallel)
#           library(CombMSC)
#           source('D:/My Dropbox/R-Projekte/mergeHMM/mleHMM.R')
#           source('D:/My Dropbox/R-Projekte/mergeHMM/merging.R')
            library(mergeHMM)
          NULL
        })
        
          parallelOutput = parLapplyLB(cl = cl, X=1:Btmp,fun=tmpMle,data=x,Glist=GlistTmp,hmm=hmm0,...)

        
        stopCluster(cl)
        
        #parallelOutput = mclapply(X=1:Btmp,FUN=tmpMle,x=x,Glist=GlistTmp,hmm=hmm0,...)

        # paraOut0 <<- parallelOutput
        
        #find out which restriction yields highest likelihood
        mllkTmp = numeric(Btmp) + Inf

        for(i in 1:Btmp) try(mllkTmp[i] <- parallelOutput[[i]]$mllk)

        if(all(mllkTmp==Inf)){
          cat("All optimizations in this step failed.")
          break
        }
        
        jTmp = which.min(mllkTmp)

        #set partition which yields highest likelihood as new one
        Gpath[[l]] = GlistTmp[[jTmp]]
        
        modelPath[[l]] = parallelOutput[[jTmp]]
        
        mllkPath[l] = modelPath[[l]]$mllk

        likelihoodRatioPath[l] = 2*(modelPath[[l]]$mllk - modelPath[[m]]$mllk)
        
        rNow = length(Gpath[[l]])
        
        pValuePath[l] = pchisq(q=likelihoodRatioPath[l],
                               df=m^2-2*m-rNow^2+2*rNow,lower.tail=FALSE)
        

        #check if H_0: lambda_Gpath[[l]](gamma) = gamma can be rejected
        if(pValuePath[l]<alpha){
          finalModel = l+1
          break
        }
      }
    }
    if(details==FALSE){
      return(list(pValuePath = rev(pValuePath),
                  finalModel = modelPath[[finalModel]]$mle,
                  finalG=Gpath[[finalModel]]))
    }
    else{
      return(list(modelPath=rev(modelPath),
                  Gpath=rev(Gpath),
                  loglikePath=rev(-mllkPath),
                  likelihoodRatioPath=rev(likelihoodRatioPath),
                  pValuePath=rev(pValuePath),
                  finalModel = modelPath[[finalModel]]$mle,
                  finalG=Gpath[[finalModel]]))  
    }
  }
  
#   #for each partition the tpm is mapped with reduceTpm to a tpm
#   #fulfilling gamma = \lambda_G(gamma)
#   if(estimateRestrictedModels == FALSE){
#     
#     #use EM-Algo. fit-model as mle
#     fullModel = HMMFit(obs=x, dis="NORMAL", nStates=m, ...)
#     
#     
#     #track all models during reducing steps
#     modelPath = list()
#     modelPath[[m]] = fullModel
#     
#     
#     #track path of log-likelihoods
#     mllkPath = NA * numeric(m) 
#     mllkPath[m] = -fullModel$LLH
#     
#     
#     gamma0 = fullModel$HMM$transMat
#     
#     if(distribution=="normal"){
#       mu0 = fullModel$HMM$distribution$mean
#       sigma0 = sqrt(fullModel$HMM$distribution$var)
# 
#     } 
#     else{
#       mu0 = matrix(0,ncol=m,nrow=mDim)*NA
#       sigma0 = array(0,dim=c(mDim,mDim,m))*NA
#       for(i in 1:m){
#         mu0[,i] = fullModel$HMM$distribution$mean[[i]]
#         sigma0[,,i] = fullModel$HMM$distribution$cov[[i]]
#       }
#     }
#     
#     stateParam0 = list(mu=mu0,sigma=sigma0)
#     
#     
#     #set a function which calculates all mle's of a partition list
#     tmpReduceHmm = function(i,Glist,x,distribution){
#       modelTmp = list(stateParameters = stateParam0,gamma = NA,mllk=NA)
#       modelTmp$gamma = reduceTpm(gamma=gamma0,G=Glist[[i]])
#       modelTmp$mllk = HMMlogLike(stateParameters=stateParam0,gamma=modelTmp$gamma,
#                                  x=x,distribution=distribution,...)
#       return(modelTmp)
#     }
#     
# 
#     
#     for(l in (m-1):1){
#   
#       #last step: iid model
#       if(l==1){        
#         modelPath[[l]] = list(stateParameters = stateParam0,gamma = NA,mllk=NA)
#         modelPath[[l]]$gamma = reduceTpm(gamma=gamma0,G=list(1:m))
#         modelPath[[l]]$mllk = HMMlogLike(stateParameters=stateParam0,gamma=modelPath[[l]]$gamma,
#                                          x=x,distribution=distribution,...)
#         Gpath[[l]] = list(1:m) 
#         mllkPath[l] = modelPath[[l]]$mllk
#         
#         likelihoodRatioPath[l] = 2*(modelPath[[l]]$mllk - mllkPath[m])
#         rNow = length(Gpath[[l]])      
#         pValuePath[l] = pchisq(q=likelihoodRatioPath[l],
#                                df=m^2-2*m-rNow^2+2*rNow,lower.tail=FALSE)
#         
#         #check if H_0: lambda_Gpath[[1]](gamma) = gamma can be rejected
#         if(pValuePath[l]<alpha){
#           finalModel = l+1
#         }
#         else
#           finalModel = l        
#       }# end if l==1
#       else{
#         #calculate for all one-level-lower-partitions the mles under the restriction
#         GlistTmp = getAllSmallerPartitions(G=Gpath[[l+1]])
#         Btmp = length(GlistTmp)
#         parallelOutput = lapply(X=1:Btmp,tmpReduceHmm,x=x,distribution=distribution,
#                                   Glist=GlistTmp)
#         
#         #print(parallelOutput)
#         
#         #find out which restriction yields highest likelihood
#         mllkTmp = numeric(Btmp)
#         for(i in 1:Btmp) mllkTmp[i] = parallelOutput[[i]]$mllk
#         jTmp = which.min(mllkTmp)
#         
#         #set partition which yields highest likelihood as new one
#         Gpath[[l]] = GlistTmp[[jTmp]]
#         
#         modelPath[[l]] = parallelOutput[[jTmp]]
#         
#         mllkPath[l] = modelPath[[l]]$mllk
#         
#         likelihoodRatioPath[l] = 2*(modelPath[[l]]$mllk - mllkPath[m])
#         
#         rNow = length(Gpath[[l]])
#         
#         pValuePath[l] = pchisq(q=likelihoodRatioPath[l],
#                                df=m^2-2*m-rNow^2+2*rNow,lower.tail=FALSE)
# 
#         #check if H_0: lambda_Gpath[[l]](gamma) = gamma can be rejected
#         if(pValuePath[l]<alpha){
#           finalModel = l+1
#           break
#         }
#       }
#     }
#     if(details==FALSE){
#       return(list(finalModel = modelPath[[finalModel]],finalG=Gpath[[finalModel]]))
#     }
#     else{
#       return(list(modelPath=modelPath,Gpath=Gpath,loglikePath=-mllkPath,
#                   likelihoodRatioPath=likelihoodRatioPath,pValuePath=pValuePath,
#                   finalModel = modelPath[[finalModel]],finalG=Gpath[[finalModel]])) 
#     }
#         
#   }
  
}
