R/AllClasses.R

In STAN: The Genomic STate ANnotation Package

#'
#' This class is a generic container for different emission functions of Hidden Markov Models.
#' 
#' @slot type The type of emission function c('Gaussian').
#' @slot parameters A list containing the the parameters for each state.
#' @slot dim Number of dimensions.
#' @slot nStates The number of states.
#' @examples
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
# 
#' HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means))
#' @exportClass HMMEmission
.HMMEmission <- setClass("HMMEmission", representation(type = "character", 
    parameters = "list", dim = "integer", nStates = "numeric"))

#' Internal function used to inititialize an HMMEmission object.
#' 
#' @keywords internal
#' @noRd
checkParameters = function(type, parameters, nStates) {
    val = FALSE
    if (class(parameters) != "list") 
        stop("Emission parameters must be submitted as a list!")
    if (type != "JointlyIndependent" & type != "null") {
        for (i in 1:length(parameters)) {
            if (!length(parameters[[i]]) == nStates & !names(parameters)[i] %in% 
                c("sharedCov", "updateCov", "reverseComplementary")) 
                stop("Number of states does not match length of parameters")
        }
    }
    
    
    if (type == "Gaussian") {
        if (!(all(c("mu", "cov") %in% names(parameters)))) 
            stop("Gaussian emissions must contain parameters \"mu\" and \"cov\"!")
        means = parameters$mu
        covs = parameters$cov
        if (!(class(means) == "list" & class(covs) == "list")) 
            stop("Means and covariances need to be of class list!")
        if (!(nStates == length(means) & nStates == length(covs))) 
            stop("Number of states does not match dimensions of means and covariances!")
        for (i in 1:length(means)) {
            stopifnot(class(covs[[i]])[1] == "matrix")
            stopifnot(dim(covs[[i]])[1] == dim(covs[[i]])[2] & length(means[[i]]) == 
                dim(covs[[i]])[1])
            stopifnot((is.numeric(covs[[i]])) & (is.numeric(means[[i]])))
            if (!all(diag(covs[[i]]) > 0)) 
                stop("Variance is negative gaussian emissions of state ", 
                  i, ".", sep = "")
        }
        val = TRUE
    } else if (type == "Bernoulli") {
        if (class(parameters$p) != "list") 
            stop("Bernoulli emission parameters must be submitted as a list (one entry per state)!")
        if (!"p" %in% names(parameters)) 
            stop("Bernoulli emission parameters must contain entry p storing emission probabilities!")
        myLens = sapply(parameters$p, function(x) length(x) == 1)
        if (!all(myLens)) 
            stop(paste("Dimensionality of Bernoulli probabilities unequal to one for states: ", 
                which(!myLens), sep = ""))
        myP = unlist(parameters$p)
        if (!all(myP <= 1 | myP >= 0)) 
            stop("Bernoulli probabilities are not 0 <= p <= 1.")
        val = TRUE
    } else if (type == "Poisson") {
        if (!"lambda" %in% names(parameters)) 
            stop("Poisson emission parameters must define lambda!")
        if (class(parameters$lambda) != "list") 
            stop("Poisson emission parameters must be submitted as a list (one entry per state)!")
        myLens = sapply(parameters$lambda, function(x) length(x) == 1)
        if (!all(myLens)) 
            stop(paste("Dimensionality of Poisson probabilities unequal to one for states: ", 
                which(!myLens), "\n", sep = ""))
        myLambda = unlist(parameters$lambda)
        # print(paste("any na in labdas:", is.na(myLambda)))
        if (!all(myLambda > 0)) 
            stop("Poisson means are not > 0. Line 83")
        
        val = TRUE
    } else if (type == "PoissonLogNormal") {
        if (!all(c("mu", "sigma") %in% names(parameters))) 
            stop("PoissonLogNormal emission parameters must define mu and sigma!")
        if (class(parameters$mu) != "list") 
            stop("PoissonLogNormal emission parameters must be submitted as a list (one entry per state)!")
        if (class(parameters$sigma) != "list") 
            stop("PoissonLogNormal emission parameters must be submitted as a list (one entry per state)!")
        myLens = sapply(parameters$mu, function(x) length(x) == length(parameters$mu[[1]]))
        if (!all(myLens)) 
            stop(paste("Unequal lengths: ", which(!myLens), "\n", sep = ""))
        mySigmas = unlist(parameters$sigma)
        if (!all(mySigmas > 0)) 
            stop("PoissonLogNormal sigmas are not > 0.")
        
        val = TRUE
    } else if (type == "NegativeBinomial") {
        if (!all(c("mu", "size", "pi") %in% names(parameters))) 
            stop("Negative Binomial emission parameters must contain entry mu, size, sizeFactor and zero inflation frequency pi!")
        if (class(parameters$mu) != "list") 
            stop("NegativeBinomial emission parameters must be submitted as a list (one entry per state)!")
        if (class(parameters$size) != "list") 
            stop("NegativeBinomial emission parameters must be submitted as a list (one entry per state)!")
        myLens = sapply(parameters$mu, function(x) length(x) == 1)
        if (!all(myLens)) 
            stop(paste("Dimensionality of NegativeBinomial probabilities unequal to one for states: ", 
                which(!myLens), "\n", sep = ""))
        myLens = sapply(parameters$size, function(x) length(x) == 1)
        # print(myLens)
        if (!all(myLens)) 
            stop(paste("Dimensionality of NegativeBinomial probabilities unequal to one for states: ", 
                which(!myLens), "\n", sep = ""))
        mySizes = unlist(parameters$size)
        if (!all(mySizes > 0)) 
            stop("NegativeBinomial size are not > 0.")
        
        val = TRUE
    } else if (type == "Multinomial") {
        if (!all(c("p") %in% names(parameters))) 
            stop("Multinomial emission parameters must contain entry p storing emission probabilities!")
        if (class(parameters$p) != "list") 
            stop("Multinomial emission parameters must be submitted as a list (one entry per state)!")
        myLens = sapply(parameters$p, function(x) length(x) == length(parameters$p[[1]]))
        if (!all(myLens)) 
            stop(paste("Unequal lengths: ", which(!myLens), sep = ""))
        if (!all(abs(sapply(parameters$p, sum) - 1) < 1e-06)) 
            stop("Multinomial probabilities must sum up to 1.")
        myP = unlist(parameters$p)
        if (!all(myP <= 1 | myP >= 0)) 
            stop("Multinomial probabilities are not 0 <= p <= 1.")
        if (!all(sapply(parameters$p, class) == "numeric")) 
            stop("Multinomial probabilities must be vectors.")
        
        val = TRUE
    } else if (type == "JointlyIndependent") {
        if (class(parameters$emissions) != "list") 
            stop("Jointly Independent emission parameters must be submitted as a list (one entry per emission function)!\n")
        myClasses = sapply(parameters$p, function(x) class(x) == "HMMEmission")
        if (!all(myClasses)) 
            stop(paste("Parameters of Jointly Independent emissions are not of class HMMEmission: ", 
                which(!myLens), "\n", sep = ""))
        if (unique(sapply(parameters$emissions, function(x) x@nStates)) != 
            nStates) 
            stop("Jointly Independent emission must have same state number as individual emission!\n")
        
        val = TRUE
    } else if (type == "null") {
        # cat('You are not using an Emission Function! Must supply fixed
        # Emission probabilities!')
        val = TRUE
    } else {
        stop("Unknown emission function specified: ", type, "\n", sep = "")
    }
    val
}

#' Internal function used to convert Bernoulli representation in Independent Model
#' 
#' @keywords internal
#' @noRd
changeParams = function(type, parameters) {
    if (type == "JointlyIndependent") {
        myDs = sapply(parameters$emissions, function(x) x@dim)
        if (!all(myDs > 0)) 
            stop("Dimensions must be greater than 0!")
        myDs = cumsum(myDs)
        emissionDims = list()
        for (i in 1:length(myDs)) {
            if (i == 1) {
                emissionDims[[i]] = 1:myDs[i]
            } else {
                emissionDims[[i]] = (myDs[i - 1] + 1):myDs[i]
            }
            
        }
        parameters$emissionDim = emissionDims
        parameters$types = sapply(parameters$emissions, function(x) x@type)
    } else if (type == "Multinomial") {
        if (!"reverseComplementary" %in% names(parameters)) {
            parameters$reverseComplementary = as.integer(1:getDim(type, 
                parameters))
            # print('fuck it')
        } else {
            parameters$reverseComplementary = as.integer(parameters$reverseComplementary)
        }
    }
    return(parameters)
    
}



#' Get dimensions of the emission function
#' 
#' @keywords internal
#' @noRd
getDim = function(type, parameters) {
    dim = NULL
    if (type == "Gaussian") {
        dim = as.integer(length(parameters$mu[[1]]))
    }
    if (type == "Bernoulli") {
        dim = as.integer(length(parameters$p[[1]]))
    }
    if (type == "Poisson") {
        dim = as.integer(length(parameters$lambda[[1]]))
    }
    if (type == "Multinomial") {
        dim = as.integer(length(parameters$p[[1]]))
    }
    if (type == "NegativeBinomial") {
        dim = as.integer(length(parameters$mu[[1]]))
    }
    if (type == "PoissonLogNormal") {
        dim = as.integer(length(parameters$mu[[1]]))
    }
    if (type == "JointlyIndependent") {
        dim = as.integer(sum(sapply(parameters$emissions, function(x) x@dim)))
    }
    if (type == "null") {
        if (is.null(parameters$dim)) 
            stop("Must supply number of dimensions in parameter list if emission type is 'null'!")
        dim = parameters$dim
    }
    dim
}

#'
#' This function creates a HMMEmission object.
#' 
#' @title Create a HMMEmission object
#' 
#' @param type The type of emission function c('Gaussian').
#' @param parameters A list containing the the parameters for each state.
#' @param nStates The number of states.
#' 
#' @return HMMEmission
#' @examples
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
# 
#' HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means))
#' @export HMMEmission
HMMEmission = function(type = character(), parameters = list(), nStates = numeric()) {
    if (length(type) == 0) 
        stop("Must specify type of emission distribution!")
    if (length(parameters) == 0) 
        stop("Must specify parameters of emission!")
    if (length(nStates) == 0) 
        stop("Must specify number of states for emission!")
    if (!checkParameters(type, parameters, nStates)) 
        stop("Parameters not correctly specified for distribution ", type)
    
    param = changeParams(type, parameters)
    dim = getDim(type, parameters)
    .HMMEmission(type = type, parameters = param, dim = dim, nStates = as.integer(nStates))
}


#' Prints description of HMMEmission object
#' 
#' @keywords internal
#' @noRd
setMethod(f = "show", signature = c("HMMEmission"), function(object) {
    cat("An object of class \"", class(object), "\"\n", " type: ", object@type, 
        "\n dim: ", paste(rbind(names(object@parameters$order), if (length(object@parameters$order) > 
            0) {
            ": "
        }, object@dim, if (length(object@parameters$order) > 0) {
            "; "
        })), "\n nStates: ", object@nStates, "\n", sep = "")
})


#'
#' This class is a generic container for Hidden Markov Models.
#' 
#' @slot initProb Initial state probabilities.
#' @slot transMat Transition probabilities
#' @slot emission Emission parameters as an HMMEmission object.
#' @slot nStates Number of states.
#' @slot status of the HMM. On of c('initial', 'EM').
#' @slot stateNames State names.
#' @slot dimNames Names of data tracks.
#' @slot LogLik Log likelihood of a fitted HMM.
#' 
#' @seealso \code{\linkS4class{HMMEmission}}
#' @examples 
#' 
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
#' 
#' HMM(initProb=initProb, transMat=transMat, emission=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)), nStates=nStates, status='initial')


#' @exportClass HMM
.HMM <- setClass("HMM", representation(initProb = "numeric", transMat = "matrix", 
    emission = c("list"), nStates = "numeric", status = "character", stateNames = "character", 
    dimNames = "character", LogLik = "numeric"))


#' Checks validity of bdHMM parameters
#' 
#' @keywords internal
#' @noRd
checkHMMParams = function(initProb, transMat, emission, nStates) {
    if (length(initProb) == 0) 
        stop("Must specify initial state probabilities!")
    if (nrow(transMat) == 0) 
        stop("Must specify transition matrix!")
    if (length(emission) == 0) 
        stop("Must specify emission!")
    if (length(nStates) == 0) 
        stop("Must specify number of states!")
    
    if (class(emission) == "HMMEmission") {
        emission = list(emission)
    } else if (class(emission) == "list") {
        if (!all(sapply(emission, class) == "HMMEmission")) 
            stop("Invalid object suplied in emission list. List of emissions must only contain objects of class HMMEmission!")
    } else {
        stop("Invalid object suplied in emission list.")
    }
    
    if (!all(sapply(emission, function(x) x@nStates) == nStates)) 
        stop("Number of states differ between HMMEmission and HMM object!")
    
    list(initProb = initProb, transMat = transMat, emission = emission)
}


#'
#' This function creates a HMM object.
#' @title Create a HMM object
#' 
#' @param initProb Initial state probabilities.
#' @param transMat Transition probabilities
#' @param emission Emission parameters as an HMMEmission object.
#' @param nStates Number of states.
#' @param status of the HMM. On of c('initial', 'EM').
#' @param stateNames State names.
#' @param dimNames Names of data tracks.
#' @param LogLik Log likelihood of a fitted HMM.
#' 
#' @seealso \code{\linkS4class{HMMEmission}}
#' 
#' @return HMM
#' @examples 
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
# 
#' HMM(initProb=initProb, transMat=transMat, emission=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)), nStates=nStates, status='initial')
#' 
#' @export HMM
HMM <- function(initProb = numeric(), transMat = matrix(numeric(), ncol = 1, 
    nrow = 1), emission, nStates = numeric(), status = character(), stateNames = character(), 
    dimNames = character(), LogLik = numeric()) {
    if (length(status) == 0) 
        status = "Object is an initial estimate"
    
    params = checkHMMParams(initProb, transMat, emission, nStates)
    # if(length(stateNames) == 0) stateNames = character()
    # if(length(dimNames) == 0) dimNames = NULL
    .HMM(initProb = params$initProb, transMat = params$transMat, emission = params$emission, 
        nStates = nStates, status = status, stateNames = stateNames, dimNames = dimNames, 
        LogLik = as.numeric(NA))
}




#' Prints description of HMM object
#' 
#' @keywords internal
#' @noRd
setMethod(f = "show", signature = c("HMM"), function(object) {
    cat("An object of class \"", class(object), "\" \n", sep = "")
    cat(" nStates: ", object@nStates, "\n", sep = "")
    cat(" stateNames: ", paste(object@stateNames, collapse = " "), "\n", 
        sep = "")
    cat(" dimNames: ", paste(object@dimNames, collapse = " "), "\n", sep = "")
    cat(" type of emission function(s): ", paste(unique(sapply(object@emission, 
        function(x) x@type)), collapse = ", "), "\n", sep = "")
    cat(" dim: ", paste(sum(sapply(object@emission, function(x) x@dim)), 
        collapse = ", "), "\n", sep = "")
    cat(" status: ", object@status, "\n", sep = "")
    if (!any(is.na(object@LogLik))) {
        cat(" Log-Likelihood: ", object@LogLik[length(object@LogLik)], 
            "\n", sep = "")
    }
    
})


#'
#' This class is a generic container for bidirectional Hidden Markov Models.
#' 
#' @slot initProb Initial state probabilities.
#' @slot transMat Transition probabilities
#' @slot emission Emission parameters as an HMMEmission object.
#' @slot nStates Number of states.
#' @slot status of the HMM. On of c('initial', 'EM').
#' @slot stateNames State names.
#' @slot dimNames Names of data tracks.
#' @slot LogLik Log likelihood of a fitted HMM.
#' @slot transitionsOptim There are three methods to choose from for fitting the transitions. Bidirectional transition matrices (invariant under reversal of time and direction) can be fitted using c('rsolnp', 'ipopt'). 'None' uses standard update formulas and the resulting matrix is not constrained to be bidirectional.
#' @slot directedObs An integer indicating which dimensions are directed. Undirected dimensions are 0. Directed observations must be marked as unique integer pairs. For instance c(0,0,0,0,0,1,1,2,2,3,3) contains 5 undirected observations, and thre pairs (one for each direction) of directed observations.
#' @slot dirScore Directionlity score of states of a fitted bdHMM.
#' 
#' @examples 
#' nStates = 5
#' stateNames = c('F1', 'F2', 'R1', 'R2', 'U1')
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
#' myEmission = list(d1=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)))
#' 
#' bdhmm = bdHMM(initProb=initProb, transMat=transMat, emission=myEmission, nStates=nStates, status='initial', stateNames=stateNames, transitionsOptim='none', directedObs=as.integer(0))
#' 
#' @seealso \code{\linkS4class{HMMEmission}}
#' @exportClass bdHMM
.bdHMM <- setClass("bdHMM", contains = "HMM", representation(transitionsOptim = "character", 
    directedObs = "numeric", dirScore = "numeric"))


validHMMEmission = function(myEmission) {
    checkParameters(myEmission@type, myEmission@parameters, myEmission@nStates)
}

#' Checks equality of emission of twin states
#' 
#' @keywords internal
#' @noRd
checkEmissionTwins = function(mybdHMM) {
    val = TRUE
    bdHMM.settings = list()
    bdHMM.settings$state2flag = getStateFlags(mybdHMM@stateNames)
    bdHMM.settings$rev.operation = getRevOperation(mybdHMM)
    myStateNames = mybdHMM@stateNames
    myF = myStateNames[grep("F", mybdHMM@stateNames)]
    myR = gsub("F", "R", myF)
    
    revOpEmissions = NULL
    if (mybdHMM@emission[[1]]@type == "JointlyIndependent") {
        revOpEmissions = emissionRevOp(mybdHMM@emission, type = "JointlyIndependent", 
            bdHMM.settings = bdHMM.settings, nStates = mybdHMM@nStates, 
            directedObs = mybdHMM@directedObs, dimNames = NULL, stateNames = mybdHMM@stateNames)$emissions
    } else {
        revOpEmissions = emissionRevOp(list(emissions = mybdHMM@emission), 
            type = "JointlyIndependent", bdHMM.settings = bdHMM.settings, 
            nStates = mybdHMM@nStates, directedObs = mybdHMM@directedObs, 
            dimNames = NULL, stateNames = mybdHMM@stateNames)$emissions
    }
    # print(EmissionParams(mybdHMM))
    mybdHMM@emission = revOpEmissions
    myParams = EmissionParams(mybdHMM)
    # print(myParams)
    for (p in 1:length(myParams)) {
        if (class(myParams[[p]])[1] == "list") {
            # here we have covs
            if (length(myF) > 0) {
                for (state in 1:length(myF)) {
                  if (any(myParams[[p]][[myF[state]]] != myParams[[p]][[myR[state]]])) 
                    stop("Emission parameters of twin states do not match.")
                }
            }
        } else if (class(myParams[[p]])[1] == "matrix") {
            if (length(myF) > 0) {
                for (state in 1:length(myF)) {
                  # print(myParams[[p]][1,])
                  if (any(myParams[[p]][myF[state], ] != myParams[[p]][myR[state], 
                    ])) 
                    stop("Emission parameters of twin states do not match.")
                }
                
            }
        } else {
            stop("Unknown format of parameter ", names(myParams)[p], " in output of EmissionParams().", 
                sep = "")
        }
    }
    val
}


#' Creates rev.operation vector from directedObs vector in bdHMM object
#' 
#' @keywords internal
#' @noRd
getRevOperation = function(mybdHMM) {
    rev.operation = 1:sum(sapply(mybdHMM@emission, function(x) x@dim))
    if (!all(table(mybdHMM@directedObs[mybdHMM@directedObs != 0]) == 2)) 
        stop("Matching of directed observations is not possible!")
    for (i in 1:length(mybdHMM@directedObs)) {
        if (mybdHMM@directedObs[i] == 0) {
            rev.operation[i] = i
        } else {
            myPair = which(mybdHMM@directedObs == mybdHMM@directedObs[i])
            myPair = myPair[myPair != i]
            rev.operation[i] = myPair
        }
    }
    rev.operation
}

#' Checks validity of a bdHMM object
#' 
#' @keywords internal
#' @noRd
validbdHMM = function(mybdHMM, probTol = 1) {
    myHMMBool = validHMM(mybdHMM)
    myF = grep("F", mybdHMM@stateNames)
    myR = grep("R", mybdHMM@stateNames)
    F2R = gsub("F", "R", mybdHMM@stateNames[myF])
    if (!all(F2R %in% mybdHMM@stateNames[myR])) 
        stop("No matching between twin states possible!")
    couples = rep(-1, length(mybdHMM@stateNames))
    if (length(myF) > 0) {
        for (i in 1:length(myF)) {
            myTwin = which(mybdHMM@stateNames == F2R[i])
            couples[myF[i]] = myTwin
            couples[myTwin] = myF[i]
        }
    }
    
    # check if rev.operation works
    rev.operation = getRevOperation(mybdHMM)
    # check if emission twin symmetry is satisfied
    myTwinsBool = checkEmissionTwins(mybdHMM)
    
    myReord = order(mybdHMM@stateNames)
    constraints = getConstraints(mybdHMM@stateNames[myReord])
    transMat = mybdHMM@transMat[myReord, myReord]
    initProb = mybdHMM@initProb[myReord]
    z = c()
    for (k in 1:length(constraints)) {
        first_prob = transMat[constraints[[k]][[1]][1], constraints[[k]][[1]][2]]
        sec_prob = transMat[constraints[[k]][[2]][1], constraints[[k]][[2]][2]]
        mult = initProb[constraints[[k]][[2]][1]]/initProb[constraints[[k]][[2]][2]]
        z[length(z) + 1] = first_prob - sec_prob * mult
    }
    if (any(abs(z) > probTol)) 
        warning("Transition symmetry is violated by ", max(abs(z)), ".", 
            sep = "")
    initDiff = abs(initProb %*% transMat - initProb)
    if (any(initDiff > probTol)) 
        warning("Initial state probabilities differ from steady state by ", 
            max(initDiff), ".", sep = "")
    
    TRUE
}

#' Checks validity of a HMM object
#' 
#' @keywords internal
#' @noRd
validHMM = function(myHMM) {
    nStates = myHMM@nStates
    # nStates matches between HMM and Emission
    if (!all(nStates == unique(sapply(myHMM@emission, function(x) x@nStates)))) 
        stop("Number of states do not match between HMM and emissions.")
    if (nrow(myHMM@transMat) != ncol(myHMM@transMat)) 
        stop("Transition matrix is not quadratic.")
    if (!all(dim(myHMM@transMat) == nStates)) 
        stop("Number of states does not match dimensions of transition matrix.")
    if (!length(myHMM@initProb) == nStates) 
        stop("Number of states does not match dimensions of transition matrix.")
    TRUE
}



#'
#' This function creates a bdHMM function.
#' 
#' @title Create a bdHMM object
#' 
#' @param initProb Initial state probabilities.
#' @param transMat Transition probabilities
#' @param emission Emission parameters as an HMMEmission object.
#' @param nStates Number of states.
#' @param status Status of the bdHMM. 'Initial' means that the model was not fitted yet. 'EM' means that the model was optimized using Expectation maximization.
#' @param stateNames Indicates directinality of states. States can be forward (F1, F2, ..., Fn), reverse (R1, R2, ..., Rn) or undirectional (U1, U2, ..., Um). Number of F and R states must be equal and twin states are indicated by integers in id (e.g. F1 and R1 and twins).
#' @param dimNames Names of data tracks.
#' @param transitionsOptim There are three methods to choose from for fitting the transitions. Bidirectional transition matrices (invariant under reversal of time and direction) can be fitted using c('rsolnp', 'analytical'). 'None' uses standard update formulas and the resulting matrix is not constrained to be bidirectional.
#' @param directedObs An integer indicating which dimensions are directed. Undirected dimensions are 0. Directed observations must be marked as unique integer pairs. For instance c(0,0,0,0,0,1,1,2,2,3,3) contains 5 undirected observations, and thre pairs (one for each direction) of directed observations.
#' @param dirScore Directionlity score of states of a fitted bdHMM.
#' 
#' @return bdHMM
#' @examples 
#' nStates = 5
#' stateNames = c('F1', 'F2', 'R1', 'R2', 'U1')
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
#' myEmission = list(d1=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)))
#' 
#' bdhmm = bdHMM(initProb=initProb, transMat=transMat, emission=myEmission, nStates=nStates, status='initial', stateNames=stateNames, transitionsOptim='none', directedObs=as.integer(0))
#' 
#' @seealso \code{\linkS4class{HMMEmission}}
#' @export bdHMM
bdHMM = function(initProb = numeric(), transMat = matrix(numeric(), ncol = 0, 
    nrow = 0), emission, nStates = numeric(), status = character(), stateNames = character(), 
    dimNames = character(), transitionsOptim = "analytical", directedObs = integer(), 
    dirScore = numeric()) {
    if (length(status) == 0) 
        status = "Object is an initial estimate"
    if (length(stateNames) == 0) 
        stop("Must specify stateLables for bdHMM!")
    if (length(directedObs) == 0) 
        stop("Must specify directedObs for bdHMM!")
    if (any(table(directedObs[directedObs != 0]) != 2)) 
        stop("Pairing of directed observations is not possible!")
    if (!class(emission) == "list") 
        stop("Object emission is not a list.")
    if (length(directedObs) != sum(sapply(emission, function(x) x@dim))) 
        stop("Length of directedObs does not match dimensionality of emission.")
    # print(dimNames)
    if ((length(dimNames) > 0) & length(dimNames) != sum(sapply(emission, 
        function(x) x@dim))) 
        stop("Length of dimNames does not match dimensionality of emission.")
    if (length(stateNames) != nStates) 
        stop("Length of stateNames does not match number of states.")
    
    params = checkHMMParams(initProb, transMat, emission, nStates)
    rev.operation = c()
    for (i in 1:length(directedObs)) {
        if (directedObs[i] == 0) {
            rev.operation[i] = i
        } else {
            myPair = which(directedObs == directedObs[i])
            myPair = myPair[myPair != i]
            rev.operation[i] = myPair
        }
    }
    bdHMM.settings = list()
    bdHMM.settings$state2flag = getStateFlags(stateNames)
    bdHMM.settings$rev.operation = rev.operation
    
    emissionRevOp(list(emissions = emission), type = "JointlyIndependent", 
        bdHMM.settings = bdHMM.settings, nStates = nStates, directedObs = directedObs, 
        dimNames = NULL, stateNames = stateNames)
    
    mybdHMM = .bdHMM(HMM(initProb = params$initProb, transMat = params$transMat, 
        emission = params$emission, nStates = nStates, status = status, 
        stateNames = stateNames, dimNames = dimNames, LogLik = as.numeric(NA)), 
        transitionsOptim = transitionsOptim, directedObs = directedObs, 
        dirScore = dirScore)
    validbdHMM(mybdHMM)
    
    mybdHMM
}

#' Prints description of a bdHMM
#' 
#' @keywords internal
#' @noRd
setMethod(f = "show", signature = c("bdHMM"), function(object) {
    cat("An object of class \"", class(object), "\" \n", sep = "")
    cat(" nStates: ", object@nStates, "\n", sep = "")
    cat(" stateNames: ", paste(object@stateNames, collapse = " "), "\n", 
        sep = "")
    cat(" dimNames: ", paste(object@dimNames, collapse = " "), "\n", sep = "")
    cat(" type of emission function(s): ", paste(unique(sapply(object@emission, 
        function(x) x@type)), collapse = ", "), "\n", sep = "")
    cat(" dim: ", paste(sum(sapply(object@emission, function(x) x@dim)), 
        collapse = ", "), "\n", sep = "")
    cat(" directedObs: ", paste(object@directedObs, collapse = " "), "\n", 
        sep = "")
    cat(" ", object@status, "\n", sep = "")
    # cat(' transitionsOptim: ', object@transitionsOptim, '\n', sep='')
    if (!any(is.na(object@LogLik))) {
        cat(" Log-Likelihood: ", object@LogLik[length(object@LogLik)], 
            "\n", sep = "")
    }
})





#'
#' This function returns the transition matrix of a (bd)HMM.
#' 
#' @title Get transitions of a (bd)HMM
#' 
#' @param hmm An object of class HMM or bdHMM.
# 
#' @return The transitions as a nStates x nStates matrix.
#' @usage  Transitions(hmm)
#' @seealso \code{\linkS4class{HMM}, \linkS4class{bdHMM}}
#' 
#' @examples 
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
# 
#' hmm = HMM(initProb=initProb, transMat=transMat, emission=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)), nStates=nStates, status='initial')
#' Transitions(hmm)
#' 
#' @export Transitions
Transitions = function(hmm) {
    if (!class(hmm) %in% c("HMM", "bdHMM")) 
        stop("Object must be of class HMM or bdHMM!\n")
    myTransMat = hmm@transMat
    colnames(myTransMat) = rownames(myTransMat) = StateNames(hmm)
    myTransMat
}

#'
#' This function returns the initial state probabilities of a (bd)HMM.
#' 
#' @title Get initial state probabilities of a (bd)HMM
#' 
#' @param hmm An object of class HMM or bdHMM.
# 
#' @return The initial state probabilities as a numeric vector.
#' @usage InitProb(hmm)
#' @seealso \code{\linkS4class{HMM}, \linkS4class{bdHMM}}
#' 
#' @examples 
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
# 
#' hmm = HMM(initProb=initProb, transMat=transMat, emission=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)), nStates=nStates, status='initial')
#' InitProb(hmm)
#' 
#' @export InitProb
InitProb = function(hmm) {
    if (!class(hmm) %in% c("HMM", "bdHMM")) 
        stop("Object must be of class HMM or bdHMM!\n")
    out = hmm@initProb
    names(out) = StateNames(hmm)
    out
}

#'
#' This function returns the Emission functions of a (bd)HMM.
#' 
#' @title Get Emission functions of a (bd)HMM
#' 
#' @param hmm An object of class HMM or bdHMM.
# 
#' @return An object of class HMMEmission
#' @usage Emission(hmm)
#' @seealso \code{\linkS4class{HMMEmission}}
#' 
#' @examples 
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
# 
#' hmm = HMM(initProb=initProb, transMat=transMat, emission=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)), nStates=nStates, status='initial')
#' Emission(hmm)
#' 
#' @export Emission
Emission = function(hmm) {
    if (!class(hmm) %in% c("HMM", "bdHMM")) 
        stop("Object must be of class HMM or bdHMM!\n")
    hmm@emission
}


#'
#' This function returns the parameters of emission functions of a (bd)HMM object.
#' 
#' @title Get Emission parameters of a (bd)HMM.
#' 
#' @param hmm An object of class (bd)HMM.
# 
#' @return A list containing the parameters of the Emission functions.
#' @usage EmissionParams(hmm)
#' @seealso \code{\linkS4class{HMMEmission}, \linkS4class{HMM}, \linkS4class{bdHMM}}
#' 
#' @examples 
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
# 
#' hmm = HMM(initProb=initProb, transMat=transMat, emission=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)), nStates=nStates, status='initial')
#' EmissionParams(hmm)
#' 
#' @export EmissionParams
EmissionParams = function(hmm) {
    if (!class(hmm) %in% c("HMM", "bdHMM")) 
        stop("Object must be of class HMM or bdHMM!\n")
    # lapply(hmm@emission, function(x) x@parameters)
    emission2par = list(Bernoulli = "p", Gaussian = c("mu", "cov"), NegativeBinomial = c("mu", 
        "size", "pi"), PoissonLogNormal = c("mu", "sigma"), Poisson = c("lambda"), 
        Multinomial = c("p"))
    allType = sapply(hmm@emission, function(x) x@type)
    allPars = unique(unlist(emission2par[allType]))
    myDimSets = list()
    Ds = sapply(hmm@emission, function(x) x@dim)
    myDims = cumsum(sapply(hmm@emission, function(x) x@dim))
    currOffSet = 0
    for (i in 1:length(Ds)) {
        myDimSets[[i]] = (currOffSet + 1):(currOffSet + Ds[i])
        currOffSet = currOffSet + length(myDimSets[[i]])
    }
    
    mySummary = list()
    stateNames = StateNames(hmm)
    dimName = DimNames(hmm)
    
    for (i in 1:length(hmm@emission)) {
        for (currPar in emission2par[[allType[i]]]) {
            if (currPar != "cov") {
                mySummary[[currPar]] = cbind(mySummary[[currPar]], do.call("rbind", 
                  hmm@emission[[i]]@parameters[[currPar]]))
            } else if (currPar == "cov" & all(sapply(hmm@emission[[i]]@parameters[[currPar]], 
                dim) == 1)) {
                mySummary[[currPar]] = cbind(mySummary[[currPar]], do.call("rbind", 
                  hmm@emission[[i]]@parameters[[currPar]]))
            } else {
                curr = hmm@emission[[i]]@parameters[[currPar]]
                names(curr) = stateNames
                mySummary[[currPar]] = c(mySummary[[currPar]], curr)
            }
        }
    }
    
    for (i in 1:length(hmm@emission)) {
        for (currPar in emission2par[[allType[i]]]) {
            if (currPar != "cov") {
                rownames(mySummary[[currPar]]) = stateNames
                colnames(mySummary[[currPar]]) = dimName[unlist(myDimSets[which(allType == 
                  hmm@emission[[i]]@type)])]
            } else if (currPar == "cov" & all(sapply(hmm@emission[[i]]@parameters[[currPar]], 
                dim) == 1)) {
                rownames(mySummary[[currPar]]) = stateNames
                colnames(mySummary[[currPar]]) = dimName[unlist(myDimSets[which(allType == 
                  hmm@emission[[i]]@type)])]
            } else {
                for (k in 1:length(mySummary[[currPar]])) {
                  colnames(mySummary[[currPar]][[k]]) = rownames(mySummary[[currPar]][[k]]) = dimName[unlist(myDimSets[which(allType == 
                    hmm@emission[[i]]@type)])]
                }
            }
        }
    }
    
    mySummary
}

#'
#' This function returns the names of dimensions (data tracks).
#' 
#' @title Get dimNames of a (bd)HMM
#' 
#' @param hmm An object of class HMM or bdHMM.
# 
#' @return A character vector
#' @usage DimNames(hmm)
#' 
#' @examples 
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
# 
#' hmm = HMM(dimNames="1", initProb=initProb, transMat=transMat, emission=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)), nStates=nStates, status='initial')
#' DimNames(hmm)
#' 
#' @export DimNames
DimNames = function(hmm) {
    if (!class(hmm) %in% c("HMM", "bdHMM")) 
        stop("Object must be of class HMM or bdHMM!\n")
    out = hmm@dimNames
    if (length(out) == 0) 
        out = NULL
    out
}

#'
#' This function returns the names of states.
#' 
#' @title Get stateNames of a (bd)HMM
#' 
#' @param hmm An object of class HMM or bdHMM.
# 
#' @return A character vector
#' @usage StateNames(hmm)
#' 
#' @examples 
#' nStates = 5
#' means = list(4,11,4,11,-1)
#' Sigma = lapply(list(4,4,4,4,4), as.matrix)
#' transMat = matrix(1/nStates, nrow=nStates, ncol=nStates)
#' initProb = rep(1/nStates, nStates)
# 
#' hmm = HMM(stateNames=as.character(1:5), initProb=initProb, transMat=transMat, emission=HMMEmission(type='Gaussian', parameters=list(mu=means, cov=Sigma), nStates=length(means)), nStates=nStates, status='initial')
#' StateNames(hmm)
#' 
#' @export StateNames
StateNames = function(hmm) {
    if (!class(hmm) %in% c("HMM", "bdHMM")) 
        stop("Object must be of class HMM or bdHMM!\n")
    out = hmm@stateNames
    if (length(out) == 0) 
        out = NULL
    out
}

#'
#' This function returns the Log-Likelihood of a (bd)HMM.
#' 
#' @title Get stateNames of a (bd)HMM
#' 
#' @param hmm An object of class HMM or bdHMM.
# 
#' @return Log likelihood during model fitting.
#' @usage LogLik(hmm)
#' 
#' @examples 
#' data(example)
#' hmm_ex = initHMM(observations, nStates=3, method="Gaussian")
#' hmm_fitted = fitHMM(observations, hmm_ex)
#' LogLik(hmm_fitted)
#' 
#' @export LogLik
LogLik = function(hmm) {
    if (!class(hmm) %in% c("HMM", "bdHMM")) 
        stop("Object must be of class HMM or bdHMM!\n")
    out = hmm@LogLik
    out
}

#'
#' This function returns the directionality score of a bdHMM.
#' 
#' @title Get directionality score of a bdHMM
#' 
#' @param bdhmm An object of class bdHMM.
# 
#' @return Directionality score of the bdHMM after model fitting.
#' @usage DirScore(bdhmm)
#' 
#' @examples 
#' data(example)
#' bdhmm_ex = initBdHMM(observations, dStates=3, method="Gaussian", directedObs=0)
#' 
#' # without flags
#' bdhmm_fitted_noFlags = fitHMM(observations, bdhmm_ex)
#' DirScore(bdhmm_fitted_noFlags)
#' 
#' # with flags
#' bdhmm_fitted_flags = fitHMM(observations, bdhmm_ex, dirFlags=flags)
#' DirScore(bdhmm_fitted_flags)
#' 
#' @export DirScore
DirScore = function(bdhmm) {
    if (!class(bdhmm) == "bdHMM") 
        stop("Object must be of class bdHMM!\n")
    out = bdhmm@dirScore
    if (length(out) > 0) {
        names(out) = bdhmm@stateNames
    }
    out
}


#' extract parts of HMM
#'
#' @param x A hidden Markov model.
#' @param i State ids to extract.
#' @param j Emissions to extract.
#' @param drop ...
#' @param ... ...
#' 
#' @return extract parts of HMM
setMethod("[", signature(x="HMM", i="ANY", j="ANY"), function(x, i, j, ..., drop="missing") {
    if (missing(i)) 
        i <- 1:x@nStates
    if (missing(j)) 
        j <- 1:sum(sapply(x@emission, function(y) y@dim))
    if (is.character(i)) 
        i = match(i, x@stateNames)
    if (is.character(j)) 
        j = match(j, x@dimNames)
    if (any(is.na(j))) 
        stop("Dim names do not match. HMM is not subsettable.")
    if (any(is.na(i))) 
        stop("Sate names do not match. HMM is not subsettable.")
    
    myDimSets = list()
    Ds = sapply(x@emission, function(x) x@dim)
    myDims = cumsum(sapply(x@emission, function(x) x@dim))
    currOffSet = 0
    for (k in 1:length(Ds)) {
        myDimSets[[k]] = (currOffSet + 1):(currOffSet + Ds[k])
        currOffSet = currOffSet + length(myDimSets[[k]])
    }
    if (!all(j %in% unlist(myDimSets))) 
        stop("Subscript of dimension out of bounds.")
    
    oneDim = which(sapply(myDimSets, length) == 1)
    multiDim = which(sapply(myDimSets, length) > 1)
    if (length(multiDim) > 0) {
        for (k in multiDim) {
            take = which(myDimSets[[k]] %in% j)
            for (l in 1:length(x@emission[[k]]@parameters)) {
                if (!names(x@emission[[k]]@parameters)[l] %in% c("sharedCov", 
                  "updateCov", "reverseComplementary")) {
                  for (state in 1:length(x@emission[[k]]@parameters[[l]])) {
                    if (class(x@emission[[k]]@parameters[[l]][[state]])[1] == "matrix") {
                      x@emission[[k]]@parameters[[l]][[state]] = as.matrix(x@emission[[k]]@parameters[[l]][[state]][take,take])
                    } else {
                      # print(x@emission[[k]]@parameters[[l]][[state]])
                      x@emission[[k]]@parameters[[l]][[state]] = x@emission[[k]]@parameters[[l]][[state]][take]
                    }
                  }
                }
            }
            x@emission[[k]]@dim = length(take)
        }
    }
    
    x@emission = x@emission[sapply(myDimSets, function(x) any(x %in% j))]
    
    x@stateNames = x@stateNames[i]
    for (k in 1:length(x@emission)) {
        for (l in 1:length(x@emission[[k]]@parameters)) {
            if (!names(x@emission[[k]]@parameters)[l] %in% c("sharedCov", 
                "updateCov", "reverseComplementary")) {
                x@emission[[k]]@parameters[[l]] = x@emission[[k]]@parameters[[l]][i]
            }
        }
        x@emission[[k]]@nStates = length(i)
    }
    x@initProb = x@initProb[i]
    x@transMat = x@transMat[i, i]
    x@nStates = length(i)
    for (k in 1:length(x@emission)) {
        validHMMEmission(x@emission[[k]])
    }
    validHMM(x)
    # x@dim = sum(sapply(x@emission, function(x) x@dim))
    x
})



#' extract parts of bdHMM
#'
#' @param x A bidirectional hidden Markov model.
#' @param i State ids to extract.
#' @param j Emissions to extract.
#' @param drop ...
#' @param ... ...
#' 
#' @return Extract parts of bdHMM
#' 
setMethod("[", signature(x="bdHMM", i="ANY", j="ANY"), function(x, i, j, ..., drop="missing") {
    if (missing(i)) 
        i <- 1:x@nStates
    if (missing(j)) 
        j <- 1:sum(sapply(x@emission, function(x) x@dim))
    if (length(i) != x@nStates) 
        stop("Subsetting is not allowed. Only reordering of states is possible.")
    if (length(j) != sum(sapply(x@emission, function(x) x@dim))) 
        stop("Subsetting is not allowed. Only reordering of dimensions is possible.")
    
    if (is.character(i)) 
        i = match(i, x@stateNames)
    if (is.character(j)) 
        j = match(j, x@dimNames)
    if (any(is.na(j))) 
        stop("Dim names do not match. bdHMM is not subsettable.")
    if (any(is.na(i))) 
        stop("Sate names do not match. bdHMM is not subsettable.")
    
    x = callNextMethod(x, i, j, ..., drop)
    x@directedObs = x@directedObs[j]
    x@dirScore = x@dirScore[i]
    # x@stateNames = x@stateNames[i]
    
    validbdHMM(x)
    for (k in 1:length(x@emission)) {
        validHMMEmission(x@emission[[k]])
    }
    x
})