R/muhisse.R

Defines functions print.muhisse.fit ParametersToPassMuHiSSE DownPassMuHisse GetFossilInitialsMuHiSSE GetRootProb FocalNodeProb SingleChildProb OrganizeData DevOptimizeMuHiSSE MuHiSSE

Documented in MuHiSSE

######################################################################################################################################
######################################################################################################################################
### MuHiSSE -- Expanded set of MuSSE models for examining diversification in relation to multiple trait evolution
######################################################################################################################################
######################################################################################################################################

MuHiSSE <- function(phy, data, f=c(1,1,1,1), turnover=c(1,2,3,4), eps=c(1,2,3,4), hidden.states=FALSE, trans.rate=NULL, condition.on.survival=TRUE, root.type="madfitz", root.p=NULL, includes.fossils=FALSE, k.samples=NULL, sann=TRUE, sann.its=1000, bounded.search=TRUE, max.tol=.Machine$double.eps^.50, starting.vals=NULL, turnover.upper=10000, eps.upper=3, trans.upper=100, restart.obj=NULL, ode.eps=0, dt.threads=1){
    
    ## Temporary fix for the current BUG:
    if( !is.null(phy$node.label) ) phy$node.label <- NULL
    
    if(!is.null(root.p)) {
        if(hidden.states ==TRUE){
            if( length( root.p ) == 4 ){
                root.p <- rep(root.p, 4)
                root.p <- root.p / sum(root.p)
                warning("For hidden states, you need to specify the root.p for all hidden states. We have adjusted it so that there's equal chance for among all hidden states.")
            } else{
                root.p.new <- numeric(32)
                root.p.new[1:length(root.p)] <- root.p
                root.p <- root.p.new
                root.p <- root.p / sum(root.p)
            }
        }else{
            ## All good:
            root.p <- root.p / sum(root.p)
            
        }
    }
    
    setDTthreads(threads=dt.threads)
    
    if(any(f == 0)){
        f[which(f==0)] <- 1
    }
    
    #if(!is.ultrametric(phy) & includes.fossils == FALSE){
    #    warning("Tree is not ultrametric. Used force.ultrametric() function to coerce the tree to be ultrametric - see note above.")
    #    edge_details <- GetEdgeDetails(phy, includes.intervals=FALSE, intervening.intervals=NULL)
    #    if(any(edge_details$type == "extinct_tip")){
    #        phy <- force.ultrametric(phy)
    #    }
    #}

    if(sann == FALSE & is.null(starting.vals)){
        warning("You have chosen to rely on the internal starting points that generally work but does not guarantee finding the MLE.")
    }
    
    if(!root.type == "madfitz" & !root.type == "herr_als"){
        stop("Check that you specified a proper root.type option. Options are 'madfitz' or 'herr_als'. See help for more details.", call.=FALSE)
    }
    
    if(is.null(trans.rate)){
        stop("Rate matrix needed. See TransMatMakerMuHiSSE() to create one.")
    }
    
    if(hidden.states == TRUE & dim(trans.rate)[1]<4){
        stop("You chose a hidden state but this is not reflected in the transition matrix")
    }
    
    ## Return error message if the data is not in the correct format.
    if( !inherits(data, what = c("matrix", "data.frame")) ){
        stop("'data' needs to be a matrix or data.frame with 3 columns. See help.")
    }
    if( !ncol( data ) == 3 ){
        stop("'data' needs to be a matrix or data.frame with 3 columns. See help.")
    }
    
    ## Check if 'hidden.states' parameter is congruent with the turnover vector:
    if( length(turnover) > 4 & !hidden.states ){
        stop("Turnover has more than 4 elements but 'hidden.states' was set to FALSE. Please set 'hidden.states' to TRUE if the model include more than one rate class.")
    }
    if( length(turnover) == 4 & hidden.states ){
        stop("Turnover has only 4 elements but 'hidden.states' was set to TRUE. Please set 'hidden.states' to FALSE if the model does not include hidden rate classes.")
    }
    
    pars <- numeric(384)
    
    if(dim(trans.rate)[2]==4){
        rate.cats <- 1
        pars.tmp <- turnover
        eps.tmp <- eps
        eps.tmp[which(eps.tmp > 0)] = (eps.tmp[which( eps.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, eps.tmp)
        trans.tmp <- c(trans.rate["(00)", "(01)"], trans.rate["(00)", "(10)"], trans.rate["(00)", "(11)"], trans.rate["(01)", "(00)"], trans.rate["(01)", "(10)"], trans.rate["(01)", "(11)"], trans.rate["(10)", "(00)"], trans.rate["(10)", "(01)"], trans.rate["(10)", "(11)"], trans.rate["(11)", "(00)"], trans.rate["(11)", "(01)"], trans.rate["(11)", "(10)"])
        trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
        category.rates.unique <- 0
        pars.tmp <- c(pars.tmp, trans.tmp)
        pars[1:20] <- pars.tmp
    }
    
    if(dim(trans.rate)[2]==8){
        rate.cats <- 2
        pars.tmp <- turnover
        eps.tmp <- eps
        eps.tmp[which(eps.tmp > 0)] = (eps.tmp[which( eps.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, eps.tmp)
        for.late.adjust <- max(pars.tmp)
        rows <- c("(00A)", "(00A)", "(00A)", "(01A)", "(01A)", "(01A)", "(10A)", "(10A)", "(10A)", "(11A)", "(11A)", "(11A)", "(00B)", "(00B)", "(00B)", "(01B)", "(01B)", "(01B)", "(10B)", "(10B)", "(10B)", "(11B)", "(11B)", "(11B)")
        cols <- c("(01A)", "(10A)", "(11A)", "(00A)", "(10A)", "(11A)", "(00A)", "(01A)", "(11A)", "(00A)", "(01A)", "(10A)", "(01B)", "(10B)", "(11B)", "(00B)", "(10B)", "(11B)", "(00B)", "(01B)", "(11B)", "(00B)", "(01B)", "(10B)")
        trans.tmp <- trans.rate[cbind(rows,cols)]
        trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, trans.tmp)
        rows <- c("(00A)", "(01A)", "(10A)", "(11A)", "(00B)", "(01B)", "(10B)", "(11B)")
        cols <- c("(00B)", "(01B)", "(10B)", "(11B)", "(00A)", "(01A)", "(10A)", "(11A)")
        category.tmp <- trans.rate[cbind(rows,cols)]
        category.tmp[category.tmp==0] <- NA
        category.rate.shift <- category.tmp + for.late.adjust
        category.rate.shift[is.na(category.rate.shift)] <- 0
        category.rate.shiftA <- c(category.rate.shift[1], rep(0,6), category.rate.shift[2], rep(0,6), category.rate.shift[3], rep(0,6), category.rate.shift[4], rep(0,6))
        category.rate.shiftB <- c(category.rate.shift[5], rep(0,6), category.rate.shift[6], rep(0,6), category.rate.shift[7], rep(0,6), category.rate.shift[8], rep(0,6))
        pars.tmp <- c(turnover[1:4], eps.tmp[1:4], trans.tmp[1:12], category.rate.shiftA, turnover[5:8], eps.tmp[5:8], trans.tmp[13:24], category.rate.shiftB)
        pars[1:length(pars.tmp)] <- pars.tmp
    }
    
    if(dim(trans.rate)[2]==12){
        rate.cats <- 3
        pars.tmp <- turnover
        eps.tmp <- eps
        eps.tmp[which(eps.tmp > 0)] = (eps.tmp[which( eps.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, eps.tmp)
        for.late.adjust <- max(pars.tmp)
        rows <- c("(00A)", "(00A)", "(00A)", "(01A)", "(01A)", "(01A)", "(10A)", "(10A)", "(10A)", "(11A)", "(11A)", "(11A)", "(00B)", "(00B)", "(00B)", "(01B)", "(01B)", "(01B)", "(10B)", "(10B)", "(10B)", "(11B)", "(11B)", "(11B)", "(00C)", "(00C)", "(00C)", "(01C)", "(01C)", "(01C)", "(10C)", "(10C)", "(10C)", "(11C)", "(11C)", "(11C)")
        cols <- c("(01A)", "(10A)", "(11A)", "(00A)", "(10A)", "(11A)", "(00A)", "(01A)", "(11A)", "(00A)", "(01A)", "(10A)", "(01B)", "(10B)", "(11B)", "(00B)", "(10B)", "(11B)", "(00B)", "(01B)", "(11B)", "(00B)", "(01B)", "(10B)", "(01C)", "(10C)", "(11C)", "(00C)", "(10C)", "(11C)", "(00C)", "(01C)", "(11C)", "(00C)", "(01C)", "(10C)")
        trans.tmp <- trans.rate[cbind(rows,cols)]
        trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, trans.tmp)
        rows <- c("(00A)", "(00A)", "(01A)", "(01A)", "(10A)", "(10A)", "(11A)", "(11A)", "(00B)", "(00B)", "(01B)", "(01B)", "(10B)", "(10B)", "(11B)", "(11B)", "(00C)", "(00C)", "(01C)", "(01C)", "(10C)", "(10C)", "(11C)", "(11C)")
        cols <- c("(00B)", "(00C)", "(01B)", "(01C)", "(10B)", "(10C)", "(11B)", "(11C)", "(00A)", "(00C)", "(01A)", "(01C)", "(10A)", "(10C)", "(11A)", "(11C)", "(00A)", "(00B)", "(01A)", "(01B)", "(10A)", "(10B)", "(11A)", "(11B)")
        category.tmp <- trans.rate[cbind(rows,cols)]
        category.tmp[category.tmp==0] <- NA
        category.rate.shift <- category.tmp + for.late.adjust
        category.rate.shift[is.na(category.rate.shift)] <- 0
        category.rate.shiftA <- c(category.rate.shift[1:2], rep(0,5), category.rate.shift[3:4], rep(0,5), category.rate.shift[5:6], rep(0,5), category.rate.shift[7:8], rep(0,5))
        category.rate.shiftB <- c(category.rate.shift[9:10], rep(0,5), category.rate.shift[11:12], rep(0,5), category.rate.shift[13:14], rep(0,5), category.rate.shift[15:16], rep(0,5))
        category.rate.shiftC <- c(category.rate.shift[17:18], rep(0,5), category.rate.shift[19:20], rep(0,5), category.rate.shift[21:22], rep(0,5), category.rate.shift[23:24], rep(0,5))
        pars.tmp <- c(turnover[1:4], eps.tmp[1:4], trans.tmp[1:12], category.rate.shiftA, turnover[5:8], eps.tmp[5:8], trans.tmp[13:24], category.rate.shiftB, turnover[9:12], eps.tmp[9:12], trans.tmp[25:36], category.rate.shiftC)
        pars[1:length(pars.tmp)] <- pars.tmp
    }
    
    if(dim(trans.rate)[2]==16){
        rate.cats <- 4
        pars.tmp <- turnover
        eps.tmp <- eps
        eps.tmp[which(eps.tmp > 0)] = (eps.tmp[which( eps.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, eps.tmp)
        for.late.adjust <- max(pars.tmp)
        rows <- c("(00A)", "(00A)", "(00A)", "(01A)", "(01A)", "(01A)", "(10A)", "(10A)", "(10A)", "(11A)", "(11A)", "(11A)", "(00B)", "(00B)", "(00B)", "(01B)", "(01B)", "(01B)", "(10B)", "(10B)", "(10B)", "(11B)", "(11B)", "(11B)", "(00C)", "(00C)", "(00C)", "(01C)", "(01C)", "(01C)", "(10C)", "(10C)", "(10C)", "(11C)", "(11C)", "(11C)", "(00D)", "(00D)", "(00D)", "(01D)", "(01D)", "(01D)", "(10D)", "(10D)", "(10D)", "(11D)", "(11D)", "(11D)")
        cols <- c("(01A)", "(10A)", "(11A)", "(00A)", "(10A)", "(11A)", "(00A)", "(01A)", "(11A)", "(00A)", "(01A)", "(10A)", "(01B)", "(10B)", "(11B)", "(00B)", "(10B)", "(11B)", "(00B)", "(01B)", "(11B)", "(00B)", "(01B)", "(10B)", "(01C)", "(10C)", "(11C)", "(00C)", "(10C)", "(11C)", "(00C)", "(01C)", "(11C)", "(00C)", "(01C)", "(10C)", "(01D)", "(10D)", "(11D)", "(00D)", "(10D)", "(11D)", "(00D)", "(01D)", "(11D)", "(00D)", "(01D)", "(10D)")
        trans.tmp <- trans.rate[cbind(rows,cols)]
        trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, trans.tmp)
        rows <- c("(00A)", "(00A)", "(00A)", "(01A)", "(01A)", "(01A)", "(10A)", "(10A)", "(10A)", "(11A)", "(11A)", "(11A)", "(00B)", "(00B)", "(00B)", "(01B)", "(01B)", "(01B)", "(10B)", "(10B)", "(10B)", "(11B)", "(11B)", "(11B)", "(00C)", "(00C)", "(00C)", "(01C)", "(01C)", "(01C)", "(10C)", "(10C)", "(10C)", "(11C)", "(11C)", "(11C)", "(00D)", "(00D)", "(00D)", "(01D)", "(01D)", "(01D)", "(10D)", "(10D)", "(10D)", "(11D)", "(11D)", "(11D)")
        cols <- c("(00B)", "(00C)", "(00D)", "(01B)", "(01C)", "(01D)", "(10B)", "(10C)", "(10D)", "(11B)", "(11C)", "(11D)", "(00A)", "(00C)", "(00D)", "(01A)", "(01C)", "(01D)", "(10A)", "(10C)", "(10D)", "(11A)", "(11C)", "(11D)", "(00A)", "(00B)", "(00D)", "(01A)", "(01B)", "(01D)", "(10A)", "(10B)", "(10D)", "(11A)", "(11B)", "(11D)", "(00A)", "(00B)", "(00C)", "(01A)", "(01B)", "(01C)", "(10A)", "(10B)", "(10C)", "(11A)", "(11B)", "(11C)")
        category.tmp <- trans.rate[cbind(rows,cols)]
        category.tmp[category.tmp==0] <- NA
        category.rate.shift <- category.tmp + for.late.adjust
        category.rate.shift[is.na(category.rate.shift)] <- 0
        category.rate.shiftA <- c(category.rate.shift[1:3], rep(0,4), category.rate.shift[4:6], rep(0,4), category.rate.shift[7:9], rep(0,4), category.rate.shift[10:12], rep(0,4))
        category.rate.shiftB <- c(category.rate.shift[13:15], rep(0,4), category.rate.shift[16:18], rep(0,4), category.rate.shift[19:21], rep(0,4), category.rate.shift[22:24], rep(0,4))
        category.rate.shiftC <- c(category.rate.shift[25:27], rep(0,4), category.rate.shift[28:30], rep(0,4), category.rate.shift[31:33], rep(0,4), category.rate.shift[34:36], rep(0,4))
        category.rate.shiftD <- c(category.rate.shift[37:39], rep(0,4), category.rate.shift[40:42], rep(0,4), category.rate.shift[43:45], rep(0,4), category.rate.shift[46:48], rep(0,4))
        pars.tmp <- c(turnover[1:4], eps.tmp[1:4], trans.tmp[1:12], category.rate.shiftA, turnover[5:8], eps.tmp[5:8], trans.tmp[13:24], category.rate.shiftB, turnover[9:12], eps.tmp[9:12], trans.tmp[25:36], category.rate.shiftC, turnover[13:16], eps.tmp[13:16], trans.tmp[37:48], category.rate.shiftD)
        pars[1:length(pars.tmp)] <- pars.tmp
    }
    
    if(dim(trans.rate)[2]==20){
        rate.cats <- 5
        pars.tmp <- turnover
        eps.tmp <- eps
        eps.tmp[which(eps.tmp > 0)] = (eps.tmp[which( eps.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, eps.tmp)
        for.late.adjust <- max(pars.tmp)
        rows <- c("(00A)", "(00A)", "(00A)", "(01A)", "(01A)", "(01A)", "(10A)", "(10A)", "(10A)", "(11A)", "(11A)", "(11A)", "(00B)", "(00B)", "(00B)", "(01B)", "(01B)", "(01B)", "(10B)", "(10B)", "(10B)", "(11B)", "(11B)", "(11B)", "(00C)", "(00C)", "(00C)", "(01C)", "(01C)", "(01C)", "(10C)", "(10C)", "(10C)", "(11C)", "(11C)", "(11C)", "(00D)", "(00D)", "(00D)", "(01D)", "(01D)", "(01D)", "(10D)", "(10D)", "(10D)", "(11D)", "(11D)", "(11D)","(00E)", "(00E)", "(00E)", "(01E)", "(01E)", "(01E)", "(10E)", "(10E)", "(10E)", "(11E)", "(11E)", "(11E)")
        cols <- c("(01A)", "(10A)", "(11A)", "(00A)", "(10A)", "(11A)", "(00A)", "(01A)", "(11A)", "(00A)", "(01A)", "(10A)", "(01B)", "(10B)", "(11B)", "(00B)", "(10B)", "(11B)", "(00B)", "(01B)", "(11B)", "(00B)", "(01B)", "(10B)", "(01C)", "(10C)", "(11C)", "(00C)", "(10C)", "(11C)", "(00C)", "(01C)", "(11C)", "(00C)", "(01C)", "(10C)", "(01D)", "(10D)", "(11D)", "(00D)", "(10D)", "(11D)", "(00D)", "(01D)", "(11D)", "(00D)", "(01D)", "(10D)", "(01E)", "(10E)", "(11E)", "(00E)", "(10E)", "(11E)", "(00E)", "(01E)", "(11E)", "(00E)", "(01E)", "(10E)")
        trans.tmp <- trans.rate[cbind(rows,cols)]
        trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, trans.tmp)
        rows <- c("(00A)","(00A)","(00A)","(00A)","(01A)","(01A)","(01A)","(01A)","(10A)","(10A)","(10A)","(10A)","(11A)","(11A)","(11A)","(11A)","(00B)","(00B)","(00B)","(00B)","(01B)","(01B)","(01B)","(01B)","(10B)","(10B)","(10B)","(10B)","(11B)","(11B)","(11B)","(11B)","(00C)","(00C)","(00C)","(00C)","(01C)","(01C)","(01C)","(01C)","(10C)","(10C)","(10C)","(10C)","(11C)","(11C)","(11C)","(11C)","(00D)","(00D)","(00D)","(00D)","(01D)","(01D)","(01D)","(01D)","(10D)","(10D)","(10D)","(10D)","(11D)","(11D)","(11D)","(11D)","(00E)","(00E)","(00E)","(00E)","(01E)","(01E)","(01E)","(01E)","(10E)","(10E)","(10E)","(10E)","(11E)","(11E)","(11E)","(11E)")
        cols <- c("(00B)","(00C)","(00D)","(00E)","(01B)","(01C)","(01D)","(01E)","(10B)","(10C)","(10D)","(10E)","(11B)","(11C)","(11D)","(11E)","(00A)","(00C)","(00D)","(00E)","(01A)","(01C)","(01D)","(01E)","(10A)","(10C)","(10D)","(10E)","(11A)","(11C)","(11D)","(11E)","(00A)","(00B)","(00D)","(00E)","(01A)","(01B)","(01D)","(01E)","(10A)","(10B)","(10D)","(10E)","(11A)","(11B)","(11D)","(11E)","(00A)","(00B)","(00C)","(00E)","(01A)","(01B)","(01C)","(01E)","(10A)","(10B)","(10C)","(10E)","(11A)","(11B)","(11C)","(11E)","(00A)","(00B)","(00C)","(00D)","(01A)","(01B)","(01C)","(01D)","(10A)","(10B)","(10C)","(10D)","(11A)","(11B)","(11C)","(11D)")
        category.tmp <- trans.rate[cbind(rows,cols)]
        category.tmp[category.tmp==0] <- NA
        category.rate.shift <- category.tmp + for.late.adjust
        category.rate.shift[is.na(category.rate.shift)] <- 0
        category.rate.shiftA <- c(category.rate.shift[1:4], rep(0,3), category.rate.shift[5:8], rep(0,3), category.rate.shift[9:12], rep(0,3), category.rate.shift[13:16], rep(0,3))
        category.rate.shiftB <- c(category.rate.shift[17:20], rep(0,3), category.rate.shift[21:24], rep(0,3), category.rate.shift[25:28], rep(0,3), category.rate.shift[29:32], rep(0,3))
        category.rate.shiftC <- c(category.rate.shift[33:36], rep(0,3), category.rate.shift[37:40], rep(0,3), category.rate.shift[41:44], rep(0,3), category.rate.shift[45:48], rep(0,3))
        category.rate.shiftD <- c(category.rate.shift[49:52], rep(0,3), category.rate.shift[53:56], rep(0,3), category.rate.shift[57:60], rep(0,3), category.rate.shift[61:64], rep(0,3))
        category.rate.shiftE <- c(category.rate.shift[65:68], rep(0,3), category.rate.shift[69:72], rep(0,3), category.rate.shift[73:76], rep(0,3), category.rate.shift[77:80], rep(0,3))
        pars.tmp <- c(turnover[1:4], eps.tmp[1:4], trans.tmp[1:12], category.rate.shiftA, turnover[5:8], eps.tmp[5:8], trans.tmp[13:24], category.rate.shiftB, turnover[9:12], eps.tmp[9:12], trans.tmp[25:36], category.rate.shiftC, turnover[13:16], eps.tmp[13:16], trans.tmp[37:48], category.rate.shiftD, turnover[17:20], eps.tmp[17:20], trans.tmp[49:60], category.rate.shiftE)
        pars[1:length(pars.tmp)] <- pars.tmp
    }
    
    if(dim(trans.rate)[2]==24){
        rate.cats <- 6
        pars.tmp <- turnover
        eps.tmp <- eps
        eps.tmp[which(eps.tmp > 0)] = (eps.tmp[which( eps.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, eps.tmp)
        for.late.adjust <- max(pars.tmp)
        rows <- c("(00A)", "(00A)", "(00A)", "(01A)", "(01A)", "(01A)", "(10A)", "(10A)", "(10A)", "(11A)", "(11A)", "(11A)", "(00B)", "(00B)", "(00B)", "(01B)", "(01B)", "(01B)", "(10B)", "(10B)", "(10B)", "(11B)", "(11B)", "(11B)", "(00C)", "(00C)", "(00C)", "(01C)", "(01C)", "(01C)", "(10C)", "(10C)", "(10C)", "(11C)", "(11C)", "(11C)", "(00D)", "(00D)", "(00D)", "(01D)", "(01D)", "(01D)", "(10D)", "(10D)", "(10D)", "(11D)", "(11D)", "(11D)", "(00E)", "(00E)", "(00E)", "(01E)", "(01E)", "(01E)", "(10E)", "(10E)", "(10E)", "(11E)", "(11E)", "(11E)", "(00F)", "(00F)", "(00F)", "(01F)", "(01F)", "(01F)", "(10F)", "(10F)", "(10F)", "(11F)", "(11F)", "(11F)")
        cols <- c("(01A)", "(10A)", "(11A)", "(00A)", "(10A)", "(11A)", "(00A)", "(01A)", "(11A)", "(00A)", "(01A)", "(10A)", "(01B)", "(10B)", "(11B)", "(00B)", "(10B)", "(11B)", "(00B)", "(01B)", "(11B)", "(00B)", "(01B)", "(10B)", "(01C)", "(10C)", "(11C)", "(00C)", "(10C)", "(11C)", "(00C)", "(01C)", "(11C)", "(00C)", "(01C)", "(10C)", "(01D)", "(10D)", "(11D)", "(00D)", "(10D)", "(11D)", "(00D)", "(01D)", "(11D)", "(00D)", "(01D)", "(10D)", "(01E)", "(10E)", "(11E)", "(00E)", "(10E)", "(11E)", "(00E)", "(01E)", "(11E)", "(00E)", "(01E)", "(10E)", "(01F)", "(10F)", "(11F)", "(00F)", "(10F)", "(11F)", "(00F)", "(01F)", "(11F)", "(00F)", "(01F)", "(10F)")
        trans.tmp <- trans.rate[cbind(rows,cols)]
        trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, trans.tmp)
        rows <- c("(00A)","(00A)","(00A)","(00A)","(00A)","(01A)","(01A)","(01A)","(01A)","(01A)","(10A)","(10A)","(10A)","(10A)","(10A)","(11A)","(11A)","(11A)","(11A)","(11A)","(00B)","(00B)","(00B)","(00B)","(00B)","(01B)","(01B)","(01B)","(01B)","(01B)","(10B)","(10B)","(10B)","(10B)","(10B)","(11B)","(11B)","(11B)","(11B)","(11B)","(00C)","(00C)","(00C)","(00C)","(00C)","(01C)","(01C)","(01C)","(01C)","(01C)","(10C)","(10C)","(10C)","(10C)","(10C)","(11C)","(11C)","(11C)","(11C)","(11C)","(00D)","(00D)","(00D)","(00D)","(00D)","(01D)","(01D)","(01D)","(01D)","(01D)","(10D)","(10D)","(10D)","(10D)","(10D)","(11D)","(11D)","(11D)","(11D)","(11D)","(00E)","(00E)","(00E)","(00E)","(00E)","(01E)","(01E)","(01E)","(01E)","(01E)","(10E)","(10E)","(10E)","(10E)","(10E)","(11E)","(11E)","(11E)","(11E)","(11E)","(00F)","(00F)","(00F)","(00F)","(00F)","(01F)","(01F)","(01F)","(01F)","(01F)","(10F)","(10F)","(10F)","(10F)","(10F)","(11F)","(11F)","(11F)","(11F)","(11F)")
        cols <- c("(00B)","(00C)","(00D)","(00E)","(00F)","(01B)","(01C)","(01D)","(01E)","(01F)","(10B)","(10C)","(10D)","(10E)","(10F)","(11B)","(11C)","(11D)","(11E)","(11F)","(00A)","(00C)","(00D)","(00E)","(00F)","(01A)","(01C)","(01D)","(01E)","(01F)","(10A)","(10C)","(10D)","(10E)","(10F)","(11A)","(11C)","(11D)","(11E)","(11F)","(00A)","(00B)","(00D)","(00E)","(00F)","(01A)","(01B)","(01D)","(01E)","(01F)","(10A)","(10B)","(10D)","(10E)","(10F)","(11A)","(11B)","(11D)","(11E)","(11F)","(00A)","(00B)","(00C)","(00E)","(00F)","(01A)","(01B)","(01C)","(01E)","(01F)","(10A)","(10B)","(10C)","(10E)","(10F)","(11A)","(11B)","(11C)","(11E)","(11F)","(00A)","(00B)","(00C)","(00D)","(00F)","(01A)","(01B)","(01C)","(01D)","(01F)","(10A)","(10B)","(10C)","(10D)","(10F)","(11A)","(11B)","(11C)","(11D)","(11F)","(00A)","(00B)","(00C)","(00D)","(00E)","(01A)","(01B)","(01C)","(01D)","(01E)","(10A)","(10B)","(10C)","(10D)","(10E)","(11A)","(11B)","(11C)","(11D)","(11E)")
        category.tmp <- trans.rate[cbind(rows,cols)]
        category.tmp[category.tmp==0] <- NA
        category.rate.shift <- category.tmp + for.late.adjust
        category.rate.shift[is.na(category.rate.shift)] <- 0
        category.rate.shiftA <- c(category.rate.shift[1:5], rep(0,2), category.rate.shift[6:10], rep(0,2), category.rate.shift[11:15], rep(0,2), category.rate.shift[16:20], rep(0,2))
        category.rate.shiftB <- c(category.rate.shift[21:25], rep(0,2), category.rate.shift[26:30], rep(0,2), category.rate.shift[31:35], rep(0,2), category.rate.shift[36:40], rep(0,2))
        category.rate.shiftC <- c(category.rate.shift[41:45], rep(0,2), category.rate.shift[46:50], rep(0,2), category.rate.shift[51:55], rep(0,2), category.rate.shift[56:60], rep(0,2))
        category.rate.shiftD <- c(category.rate.shift[61:65], rep(0,2), category.rate.shift[66:70], rep(0,2), category.rate.shift[71:75], rep(0,2), category.rate.shift[76:80], rep(0,2))
        category.rate.shiftE <- c(category.rate.shift[81:85], rep(0,2), category.rate.shift[86:90], rep(0,2), category.rate.shift[91:95], rep(0,2), category.rate.shift[96:100], rep(0,2))
        category.rate.shiftF <- c(category.rate.shift[101:105], rep(0,2), category.rate.shift[106:110], rep(0,2), category.rate.shift[111:115], rep(0,2), category.rate.shift[116:120], rep(0,2))
        pars.tmp <- c(turnover[1:4], eps.tmp[1:4], trans.tmp[1:12], category.rate.shiftA, turnover[5:8], eps.tmp[5:8], trans.tmp[13:24], category.rate.shiftB, turnover[9:12], eps.tmp[9:12], trans.tmp[25:36], category.rate.shiftC, turnover[13:16], eps.tmp[13:16], trans.tmp[37:48], category.rate.shiftD, turnover[17:20], eps.tmp[17:20], trans.tmp[49:60], category.rate.shiftE, turnover[21:24], eps.tmp[21:24], trans.tmp[61:72], category.rate.shiftF)
        pars[1:length(pars.tmp)] <- pars.tmp
    }
    
    if(dim(trans.rate)[2]==28){
        rate.cats <- 7
        pars.tmp <- turnover
        eps.tmp <- eps
        eps.tmp[which(eps.tmp > 0)] = (eps.tmp[which( eps.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, eps.tmp)
        for.late.adjust <- max(pars.tmp)
        rows <- c("(00A)", "(00A)", "(00A)", "(01A)", "(01A)", "(01A)", "(10A)", "(10A)", "(10A)", "(11A)", "(11A)", "(11A)", "(00B)", "(00B)", "(00B)", "(01B)", "(01B)", "(01B)", "(10B)", "(10B)", "(10B)", "(11B)", "(11B)", "(11B)", "(00C)", "(00C)", "(00C)", "(01C)", "(01C)", "(01C)", "(10C)", "(10C)", "(10C)", "(11C)", "(11C)", "(11C)", "(00D)", "(00D)", "(00D)", "(01D)", "(01D)", "(01D)", "(10D)", "(10D)", "(10D)", "(11D)", "(11D)", "(11D)", "(00E)", "(00E)", "(00E)", "(01E)", "(01E)", "(01E)", "(10E)", "(10E)", "(10E)", "(11E)", "(11E)", "(11E)", "(00F)", "(00F)", "(00F)", "(01F)", "(01F)", "(01F)", "(10F)", "(10F)", "(10F)", "(11F)", "(11F)", "(11F)", "(00G)", "(00G)", "(00G)", "(01G)", "(01G)", "(01G)", "(10G)", "(10G)", "(10G)", "(11G)", "(11G)", "(11G)")
        cols <- c("(01A)", "(10A)", "(11A)", "(00A)", "(10A)", "(11A)", "(00A)", "(01A)", "(11A)", "(00A)", "(01A)", "(10A)", "(01B)", "(10B)", "(11B)", "(00B)", "(10B)", "(11B)", "(00B)", "(01B)", "(11B)", "(00B)", "(01B)", "(10B)", "(01C)", "(10C)", "(11C)", "(00C)", "(10C)", "(11C)", "(00C)", "(01C)", "(11C)", "(00C)", "(01C)", "(10C)", "(01D)", "(10D)", "(11D)", "(00D)", "(10D)", "(11D)", "(00D)", "(01D)", "(11D)", "(00D)", "(01D)", "(10D)", "(01E)", "(10E)", "(11E)", "(00E)", "(10E)", "(11E)", "(00E)", "(01E)", "(11E)", "(00E)", "(01E)", "(10E)", "(01F)", "(10F)", "(11F)", "(00F)", "(10F)", "(11F)", "(00F)", "(01F)", "(11F)", "(00F)", "(01F)", "(10F)", "(01G)", "(10G)", "(11G)", "(00G)", "(10G)", "(11G)", "(00G)", "(01G)", "(11G)", "(00G)", "(01G)", "(10G)")
        trans.tmp <- trans.rate[cbind(rows,cols)]
        trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, trans.tmp)
        rows <- c("(00A)","(00A)","(00A)","(00A)","(00A)","(00A)","(01A)","(01A)","(01A)","(01A)","(01A)","(01A)","(10A)","(10A)","(10A)","(10A)","(10A)","(10A)","(11A)","(11A)","(11A)","(11A)","(11A)","(11A)","(00B)","(00B)","(00B)","(00B)","(00B)","(00B)","(01B)","(01B)","(01B)","(01B)","(01B)","(01B)","(10B)","(10B)","(10B)","(10B)","(10B)","(10B)","(11B)","(11B)","(11B)","(11B)","(11B)","(11B)","(00C)","(00C)","(00C)","(00C)","(00C)","(00C)","(01C)","(01C)","(01C)","(01C)","(01C)","(01C)","(10C)","(10C)","(10C)","(10C)","(10C)","(10C)","(11C)","(11C)","(11C)","(11C)","(11C)","(11C)","(00D)","(00D)","(00D)","(00D)","(00D)","(00D)","(01D)","(01D)","(01D)","(01D)","(01D)","(01D)","(10D)","(10D)","(10D)","(10D)","(10D)","(10D)","(11D)","(11D)","(11D)","(11D)","(11D)","(11D)","(00E)","(00E)","(00E)","(00E)","(00E)","(00E)","(01E)","(01E)","(01E)","(01E)","(01E)","(01E)","(10E)","(10E)","(10E)","(10E)","(10E)","(10E)","(11E)","(11E)","(11E)","(11E)","(11E)","(11E)","(00F)","(00F)","(00F)","(00F)","(00F)","(00F)","(01F)","(01F)","(01F)","(01F)","(01F)","(01F)","(10F)","(10F)","(10F)","(10F)","(10F)","(10F)","(11F)","(11F)","(11F)","(11F)","(11F)","(11F)","(00G)","(00G)","(00G)","(00G)","(00G)","(00G)","(01G)","(01G)","(01G)","(01G)","(01G)","(01G)","(10G)","(10G)","(10G)","(10G)","(10G)","(10G)","(11G)","(11G)","(11G)","(11G)","(11G)","(11G)")
        cols <- c("(00B)","(00C)","(00D)","(00E)","(00F)","(00G)","(01B)","(01C)","(01D)","(01E)","(01F)","(01G)","(10B)","(10C)","(10D)","(10E)","(10F)","(10G)","(11B)","(11C)","(11D)","(11E)","(11F)","(11G)","(00A)","(00C)","(00D)","(00E)","(00F)","(00G)","(01A)","(01C)","(01D)","(01E)","(01F)","(01G)","(10A)","(10C)","(10D)","(10E)","(10F)","(10G)","(11A)","(11C)","(11D)","(11E)","(11F)","(11G)","(00A)","(00B)","(00D)","(00E)","(00F)","(00G)","(01A)","(01B)","(01D)","(01E)","(01F)","(01G)","(10A)","(10B)","(10D)","(10E)","(10F)","(10G)","(11A)","(11B)","(11D)","(11E)","(11F)","(11G)","(00A)","(00B)","(00C)","(00E)","(00F)","(00G)","(01A)","(01B)","(01C)","(01E)","(01F)","(01G)","(10A)","(10B)","(10C)","(10E)","(10F)","(10G)","(11A)","(11B)","(11C)","(11E)","(11F)","(11G)","(00A)","(00B)","(00C)","(00D)","(00F)","(00G)","(01A)","(01B)","(01C)","(01D)","(01F)","(01G)","(10A)","(10B)","(10C)","(10D)","(10F)","(10G)","(11A)","(11B)","(11C)","(11D)","(11F)","(11G)","(00A)","(00B)","(00C)","(00D)","(00E)","(00G)","(01A)","(01B)","(01C)","(01D)","(01E)","(01G)","(10A)","(10B)","(10C)","(10D)","(10E)","(10G)","(11A)","(11B)","(11C)","(11D)","(11E)","(11G)","(00A)","(00B)","(00C)","(00D)","(00E)","(00F)","(01A)","(01B)","(01C)","(01D)","(01E)","(01F)","(10A)","(10B)","(10C)","(10D)","(10E)","(10F)","(11A)","(11B)","(11C)","(11D)","(11E)","(11F)")
        category.tmp <- trans.rate[cbind(rows,cols)]
        category.tmp[category.tmp==0] <- NA
        category.rate.shift <- category.tmp + for.late.adjust
        category.rate.shift[is.na(category.rate.shift)] <- 0
        category.rate.shiftA <- c(category.rate.shift[1:6], rep(0,1), category.rate.shift[7:12], rep(0,1), category.rate.shift[13:18], rep(0,1), category.rate.shift[19:24], rep(0,1))
        category.rate.shiftB <- c(category.rate.shift[25:30], rep(0,1), category.rate.shift[31:36], rep(0,1), category.rate.shift[37:42], rep(0,1), category.rate.shift[43:48], rep(0,1))
        category.rate.shiftC <- c(category.rate.shift[49:54], rep(0,1), category.rate.shift[55:60], rep(0,1), category.rate.shift[61:66], rep(0,1), category.rate.shift[67:72], rep(0,1))
        category.rate.shiftD <- c(category.rate.shift[73:78], rep(0,1), category.rate.shift[79:84], rep(0,1), category.rate.shift[85:90], rep(0,1), category.rate.shift[91:96], rep(0,1))
        category.rate.shiftE <- c(category.rate.shift[97:102], rep(0,1), category.rate.shift[103:108], rep(0,1), category.rate.shift[109:114], rep(0,1), category.rate.shift[115:120], rep(0,1))
        category.rate.shiftF <- c(category.rate.shift[121:126], rep(0,1), category.rate.shift[127:132], rep(0,1), category.rate.shift[133:138], rep(0,1), category.rate.shift[139:144], rep(0,1))
        category.rate.shiftG <- c(category.rate.shift[145:150], rep(0,1), category.rate.shift[151:156], rep(0,1), category.rate.shift[157:162], rep(0,1), category.rate.shift[163:168], rep(0,1))
        pars.tmp <- c(turnover[1:4], eps.tmp[1:4], trans.tmp[1:12], category.rate.shiftA, turnover[5:8], eps.tmp[5:8], trans.tmp[13:24], category.rate.shiftB, turnover[9:12], eps.tmp[9:12], trans.tmp[25:36], category.rate.shiftC, turnover[13:16], eps.tmp[13:16], trans.tmp[37:48], category.rate.shiftD, turnover[17:20], eps.tmp[17:20], trans.tmp[49:60], category.rate.shiftE, turnover[21:24], eps.tmp[21:24], trans.tmp[61:72], category.rate.shiftF, turnover[25:28], eps.tmp[25:28], trans.tmp[73:84], category.rate.shiftG)
        pars[1:length(pars.tmp)] <- pars.tmp
    }
    
    if(dim(trans.rate)[2]==32){
        rate.cats <- 8
        pars.tmp <- turnover
        eps.tmp <- eps
        eps.tmp[which(eps.tmp > 0)] = (eps.tmp[which( eps.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, eps.tmp)
        for.late.adjust <- max(pars.tmp)
        rows <- c("(00A)", "(00A)", "(00A)", "(01A)", "(01A)", "(01A)", "(10A)", "(10A)", "(10A)", "(11A)", "(11A)", "(11A)", "(00B)", "(00B)", "(00B)", "(01B)", "(01B)", "(01B)", "(10B)", "(10B)", "(10B)", "(11B)", "(11B)", "(11B)", "(00C)", "(00C)", "(00C)", "(01C)", "(01C)", "(01C)", "(10C)", "(10C)", "(10C)", "(11C)", "(11C)", "(11C)", "(00D)", "(00D)", "(00D)", "(01D)", "(01D)", "(01D)", "(10D)", "(10D)", "(10D)", "(11D)", "(11D)", "(11D)", "(00E)", "(00E)", "(00E)", "(01E)", "(01E)", "(01E)", "(10E)", "(10E)", "(10E)", "(11E)", "(11E)", "(11E)", "(00F)", "(00F)", "(00F)", "(01F)", "(01F)", "(01F)", "(10F)", "(10F)", "(10F)", "(11F)", "(11F)", "(11F)", "(00G)", "(00G)", "(00G)", "(01G)", "(01G)", "(01G)", "(10G)", "(10G)", "(10G)", "(11G)", "(11G)", "(11G)", "(00H)", "(00H)", "(00H)", "(01H)", "(01H)", "(01H)", "(10H)", "(10H)", "(10H)", "(11H)", "(11H)", "(11H)")
        cols <- c("(01A)", "(10A)", "(11A)", "(00A)", "(10A)", "(11A)", "(00A)", "(01A)", "(11A)", "(00A)", "(01A)", "(10A)", "(01B)", "(10B)", "(11B)", "(00B)", "(10B)", "(11B)", "(00B)", "(01B)", "(11B)", "(00B)", "(01B)", "(10B)", "(01C)", "(10C)", "(11C)", "(00C)", "(10C)", "(11C)", "(00C)", "(01C)", "(11C)", "(00C)", "(01C)", "(10C)", "(01D)", "(10D)", "(11D)", "(00D)", "(10D)", "(11D)", "(00D)", "(01D)", "(11D)", "(00D)", "(01D)", "(10D)", "(01E)", "(10E)", "(11E)", "(00E)", "(10E)", "(11E)", "(00E)", "(01E)", "(11E)", "(00E)", "(01E)", "(10E)", "(01F)", "(10F)", "(11F)", "(00F)", "(10F)", "(11F)", "(00F)", "(01F)", "(11F)", "(00F)", "(01F)", "(10F)", "(01G)", "(10G)", "(11G)", "(00G)", "(10G)", "(11G)", "(00G)", "(01G)", "(11G)", "(00G)", "(01G)", "(10G)", "(01H)", "(10H)", "(11H)", "(00H)", "(10H)", "(11H)", "(00H)", "(01H)", "(11H)", "(00H)", "(01H)", "(10H)")
        trans.tmp <- trans.rate[cbind(rows,cols)]
        trans.tmp[which(trans.tmp > 0)] = (trans.tmp[which(trans.tmp > 0)] + max(pars.tmp))
        pars.tmp <- c(pars.tmp, trans.tmp)
        rows <- c("(00A)","(00A)","(00A)","(00A)","(00A)","(00A)","(00A)","(01A)","(01A)","(01A)","(01A)","(01A)","(01A)","(01A)","(10A)","(10A)","(10A)","(10A)","(10A)","(10A)","(10A)","(11A)","(11A)","(11A)","(11A)","(11A)","(11A)","(11A)","(00B)","(00B)","(00B)","(00B)","(00B)","(00B)","(00B)","(01B)","(01B)","(01B)","(01B)","(01B)","(01B)","(01B)","(10B)","(10B)","(10B)","(10B)","(10B)","(10B)","(10B)","(11B)","(11B)","(11B)","(11B)","(11B)","(11B)","(11B)","(00C)","(00C)","(00C)","(00C)","(00C)","(00C)","(00C)","(01C)","(01C)","(01C)","(01C)","(01C)","(01C)","(01C)","(10C)","(10C)","(10C)","(10C)","(10C)","(10C)","(10C)","(11C)","(11C)","(11C)","(11C)","(11C)","(11C)","(11C)","(00D)","(00D)","(00D)","(00D)","(00D)","(00D)","(00D)","(01D)","(01D)","(01D)","(01D)","(01D)","(01D)","(01D)","(10D)","(10D)","(10D)","(10D)","(10D)","(10D)","(10D)","(11D)","(11D)","(11D)","(11D)","(11D)","(11D)","(11D)","(00E)","(00E)","(00E)","(00E)","(00E)","(00E)","(00E)","(01E)","(01E)","(01E)","(01E)","(01E)","(01E)","(01E)","(10E)","(10E)","(10E)","(10E)","(10E)","(10E)","(10E)","(11E)","(11E)","(11E)","(11E)","(11E)","(11E)","(11E)","(00F)","(00F)","(00F)","(00F)","(00F)","(00F)","(00F)","(01F)","(01F)","(01F)","(01F)","(01F)","(01F)","(01F)","(10F)","(10F)","(10F)","(10F)","(10F)","(10F)","(10F)","(11F)","(11F)","(11F)","(11F)","(11F)","(11F)","(11F)","(00G)","(00G)","(00G)","(00G)","(00G)","(00G)","(00G)","(01G)","(01G)","(01G)","(01G)","(01G)","(01G)","(01G)","(10G)","(10G)","(10G)","(10G)","(10G)","(10G)","(10G)","(11G)","(11G)","(11G)","(11G)","(11G)","(11G)","(11G)","(00H)","(00H)","(00H)","(00H)","(00H)","(00H)","(00H)","(01H)","(01H)","(01H)","(01H)","(01H)","(01H)","(01H)","(10H)","(10H)","(10H)","(10H)","(10H)","(10H)","(10H)","(11H)","(11H)","(11H)","(11H)","(11H)","(11H)","(11H)")
        cols <- c("(00B)","(00C)","(00D)","(00E)","(00F)","(00G)","(00H)","(01B)","(01C)","(01D)","(01E)","(01F)","(01G)","(01H)","(10B)","(10C)","(10D)","(10E)","(10F)","(10G)","(10H)","(11B)","(11C)","(11D)","(11E)","(11F)","(11G)","(11H)","(00A)","(00C)","(00D)","(00E)","(00F)","(00G)","(00H)","(01A)","(01C)","(01D)","(01E)","(01F)","(01G)","(01H)","(10A)","(10C)","(10D)","(10E)","(10F)","(10G)","(10H)","(11A)","(11C)","(11D)","(11E)","(11F)","(11G)","(11H)","(00A)","(00B)","(00D)","(00E)","(00F)","(00G)","(00H)","(01A)","(01B)","(01D)","(01E)","(01F)","(01G)","(01H)","(10A)","(10B)","(10D)","(10E)","(10F)","(10G)","(10H)","(11A)","(11B)","(11D)","(11E)","(11F)","(11G)","(11H)","(00A)","(00B)","(00C)","(00E)","(00F)","(00G)","(00H)","(01A)","(01B)","(01C)","(01E)","(01F)","(01G)","(01H)","(10A)","(10B)","(10C)","(10E)","(10F)","(10G)","(10H)","(11A)","(11B)","(11C)","(11E)","(11F)","(11G)","(11H)","(00A)","(00B)","(00C)","(00D)","(00F)","(00G)","(00H)","(01A)","(01B)","(01C)","(01D)","(01F)","(01G)","(01H)","(10A)","(10B)","(10C)","(10D)","(10F)","(10G)","(10H)","(11A)","(11B)","(11C)","(11D)","(11F)","(11G)","(11H)","(00A)","(00B)","(00C)","(00D)","(00E)","(00G)","(00H)","(01A)","(01B)","(01C)","(01D)","(01E)","(01G)","(01H)","(10A)","(10B)","(10C)","(10D)","(10E)","(10G)","(10H)","(11A)","(11B)","(11C)","(11D)","(11E)","(11G)","(11H)","(00A)","(00B)","(00C)","(00D)","(00E)","(00F)","(00H)","(01A)","(01B)","(01C)","(01D)","(01E)","(01F)","(01H)","(10A)","(10B)","(10C)","(10D)","(10E)","(10F)","(10H)","(11A)","(11B)","(11C)","(11D)","(11E)","(11F)","(11H)","(00A)","(00B)","(00C)","(00D)","(00E)","(00F)","(00G)","(01A)","(01B)","(01C)","(01D)","(01E)","(01F)","(01G)","(10A)","(10B)","(10C)","(10D)","(10E)","(10F)","(10G)","(11A)","(11B)","(11C)","(11D)","(11E)","(11F)","(11G)")
        category.tmp <- trans.rate[cbind(rows,cols)]
        category.tmp[category.tmp==0] <- NA
        category.rate.shift <- category.tmp + for.late.adjust
        category.rate.shift[is.na(category.rate.shift)] <- 0
        category.rate.shiftA <- c(category.rate.shift[1:28])
        category.rate.shiftB <- c(category.rate.shift[29:56])
        category.rate.shiftC <- c(category.rate.shift[57:84])
        category.rate.shiftD <- c(category.rate.shift[85:112])
        category.rate.shiftE <- c(category.rate.shift[113:140])
        category.rate.shiftF <- c(category.rate.shift[141:168])
        category.rate.shiftG <- c(category.rate.shift[169:196])
        category.rate.shiftH <- c(category.rate.shift[197:224])
        pars.tmp <- c(turnover[1:4], eps.tmp[1:4], trans.tmp[1:12], category.rate.shiftA, turnover[5:8], eps.tmp[5:8], trans.tmp[13:24], category.rate.shiftB, turnover[9:12], eps.tmp[9:12], trans.tmp[25:36], category.rate.shiftC, turnover[13:16], eps.tmp[13:16], trans.tmp[37:48], category.rate.shiftD, turnover[17:20], eps.tmp[17:20], trans.tmp[49:60], category.rate.shiftE, turnover[21:24], eps.tmp[21:24], trans.tmp[61:72], category.rate.shiftF, turnover[25:28], eps.tmp[25:28], trans.tmp[73:84], category.rate.shiftG, turnover[29:32], eps.tmp[29:32], trans.tmp[85:96], category.rate.shiftH)
        pars[1:length(pars.tmp)] <- pars.tmp
    }
    
    if(includes.fossils == TRUE){
        pars <- c(pars, max(pars.tmp)+1)
    }else{
        pars <- c(pars, 0)
    }
    
    np <- max(pars)
    pars[pars==0] <- np+1
    
    cat("Initializing...", "\n")
    
    # Some new prerequisites #
    if(includes.fossils == TRUE){
        if(!is.null(k.samples)){
            psi.type <- "m+k"
            split.times <- dateNodes(phy, rootAge=max(node.depth.edgelength(phy)))[-c(1:Ntip(phy))]
            #in simulation ordering of the k samples mattered:
            k.samples <- k.samples[order(as.numeric(k.samples[,3]),decreasing=FALSE),]
            phy <- AddKNodes(phy, k.samples)
            fix.type <- GetKSampleMRCA(phy, k.samples)
            data <- AddKData(data, k.samples, muhisse=TRUE)
            no.k.samples <- length(k.samples[,1])
        }else{
            psi.type <- "m_only"
            fix.type <- NULL
            split.times <- dateNodes(phy, rootAge=max(node.depth.edgelength(phy)))[-c(1:Ntip(phy))]
            no.k.samples <- 0
        }
        gen <- FindGenerations(phy)
        data.new <- data.frame(data[,2], data[,3], row.names=data[,1])
        data.new <- data.new[phy$tip.label,]
        dat.tab <- OrganizeData(data=data.new, phy=phy, f=f, hidden.states=hidden.states, includes.fossils=includes.fossils)
        #These are all inputs for generating starting values:
        fossil.taxa <- which(dat.tab$branch.type == 1)
        fossil.ages <- dat.tab$TipwardAge[which(dat.tab$branch.type == 1)]
        n.tax.starting <- Ntip(phy)-length(fossil.taxa)-no.k.samples
    }else{
        gen <- FindGenerations(phy)
        data.new <- data.frame(data[,2], data[,3], row.names=data[,1])
        data.new <- data.new[phy$tip.label,]
        dat.tab <- OrganizeData(data=data.new, phy=phy, f=f, hidden.states=hidden.states, includes.fossils=includes.fossils)
        fossil.taxa <- NULL
        fix.type <- NULL
        psi.type <- NULL
    }
    nb.tip <- Ntip(phy)
    nb.node <- phy$Nnode
    ##########################
    
    #This is used to scale starting values to account for sampling:
    if(length(f) == 4){
        freqs <- table(apply(data.new, 1, function(x) switch(paste0(x, collapse=""), "00" = 1, "01" = 2, "10" = 3, "11" = 4, "02"=1, "20"=3, "21"=2, "12"=4, "22"=4)))
        ## Fixing the structure of the freqs vector.
        freqs.vec <- rep(0, times = 4)
        freqs.vec[as.numeric(names(freqs))] <- as.numeric( freqs )
        samp.freq.tree <- Ntip(phy) / sum(freqs.vec / f)
    }else{
        if(length(f) == Ntip(phy)){
            stop("This functionality has been temporarily removed.")
            #samp.freq.tree <- Ntip(phy) / sum(table(data.new[,1]) / mean(f[as.numeric(names(freqs))+1]))
        }else{
            stop("The vector of sampling frequencies does not match the number of tips in the tree.")
        }
    }
    if(is.null(restart.obj)){
        if(sum(eps)==0){
            if(includes.fossils == TRUE){
                stop("You input a tree that contains fossils but you are assuming no extinction. Check your function call.")
            }else{
                init.pars <- starting.point.generator(phy, 4, samp.freq.tree, yule=TRUE)
            }
        }else{
            if(includes.fossils == TRUE){
                init.pars <- starting.point.generator.fossils(n.tax=n.tax.starting, k=4, samp.freq.tree, fossil.taxa=fossil.taxa, fossil.ages=fossil.ages, no.k.samples=no.k.samples, split.times=split.times)
                psi.start <- init.pars[length(init.pars)]
            }else{
                init.pars <- starting.point.generator(phy, 4, samp.freq.tree, yule=FALSE)
            }
            if(any(init.pars[5:8] == 0)){
                init.pars[5:8] = 1e-6
            }
        }
        names(init.pars) <- NULL
        
        if(is.null(starting.vals)){
            def.set.pars <- rep(c(log(init.pars[1:4]+init.pars[5:8]), log(init.pars[5:8]/init.pars[1:4]), log(init.pars[9:20]), rep(log(.01), 28)), rate.cats)
        }else{
            ## Check if 'starting.vals' has the correct format.
            if( !length(starting.vals) %in% c(3,20) ){
                stop("Wrong length of starting.vals")
            }
            if( length(starting.vals) == 20 ){
                cat("Using developer mode for starting.vals.", "\n")
                def.set.pars <- rep(c(log(starting.vals[1:4]), log(starting.vals[5:8]), log(starting.vals[9:20]), rep(log(0.01), 28)), rate.cats)
            } else{
                def.set.pars <- rep(c(log( rep(starting.vals[1],4) ), log( rep(starting.vals[2],4) ), log( rep(starting.vals[3],12) ), rep(log(0.01), 28)), rate.cats)
            }
        }
        
        if(bounded.search == TRUE){
            upper.full <- rep(c(rep(log(turnover.upper),4), rep(log(eps.upper),4), rep(log(trans.upper),12), rep(log(10), 28)), rate.cats)
        }else{
            upper.full <- rep(21,length(def.set.pars))
        }
        
        if(includes.fossils == TRUE){
            np.sequence <- 1:(np-1)
            ip <- numeric(np-1)
            upper <- numeric(np-1)
            for(i in np.sequence){
                ip[i] <- def.set.pars[which(pars == np.sequence[i])[1]]
                upper[i] <- upper.full[which(pars == np.sequence[i])[1]]
            }
            ip <- c(ip, log(psi.start))
            upper <- c(upper, log(trans.upper))
            lower <- rep(-20, length(ip))
        }else{
            np.sequence <- 1:np
            ip <- numeric(np)
            upper <- numeric(np)
            for(i in np.sequence){
                ip[i] <- def.set.pars[which(pars == np.sequence[i])[1]]
                upper[i] <- upper.full[which(pars == np.sequence[i])[1]]
            }
            lower <- rep(-20, length(ip))
        }
    }else{
        upper <- restart.obj$upper.bounds
        lower <- restart.obj$lower.bounds
        pars <- restart.obj$index.par
        ip <- numeric(length(unique(restart.obj$index.par))-1)
        for(k in 1:length(ip)){
            ip[k] <- log(restart.obj$solution[which(restart.obj$index.par==k)][1])
        }
    }
    
    if(sann == FALSE){
        if(bounded.search == TRUE){
            cat("Finished. Beginning bounded subplex routine...", "\n")
            opts <- list("algorithm" = "NLOPT_LN_SBPLX", "maxeval" = 100000, "ftol_rel" = max.tol)
            out = nloptr(x0=ip, eval_f=DevOptimizeMuHiSSE, ub=upper, lb=lower, opts=opts, pars=pars, dat.tab=dat.tab, gen=gen, hidden.states=hidden.states, nb.tip=nb.tip, nb.node=nb.node, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, f=f, ode.eps=ode.eps, fossil.taxa=fossil.taxa, fix.type=fix.type)
            sann.counts <- NULL
            solution <- numeric(length(pars))
            solution[] <- c(exp(out$solution), 0)[pars]
            loglik = -out$objective
        }else{
            cat("Finished. Beginning subplex routine...", "\n")
            out = subplex(ip, fn=DevOptimizeMuHiSSE, control=list(reltol=max.tol, parscale=rep(0.1, length(ip))), pars=pars, dat.tab=dat.tab, gen=gen, hidden.states=hidden.states, nb.tip=nb.tip, nb.node=nb.node, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, f=f, ode.eps=ode.eps, fossil.taxa=fossil.taxa, fix.type=fix.type)
            sann.counts <- NULL
            solution <- numeric(length(pars))
            solution[] <- c(exp(out$par), 0)[pars]
            loglik = -out$value
        }
    }else{
        cat("Finished. Beginning simulated annealing...", "\n")
        out.sann = GenSA(ip, fn=DevOptimizeMuHiSSE, lower=lower, upper=upper, control=list(max.call=sann.its), pars=pars, dat.tab=dat.tab, gen=gen, hidden.states=hidden.states, nb.tip=nb.tip, nb.node=nb.node, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, f=f, ode.eps=ode.eps, fossil.taxa=fossil.taxa, fix.type=fix.type)
        sann.counts <- out.sann$counts
        cat("Finished. Refining using subplex routine...", "\n")
        opts <- list("algorithm" = "NLOPT_LN_SBPLX", "maxeval" = 100000, "ftol_rel" = max.tol)
        out <- nloptr(x0=out.sann$par, eval_f=DevOptimizeMuHiSSE, ub=upper, lb=lower, opts=opts, pars=pars, dat.tab=dat.tab, gen=gen, hidden.states=hidden.states, nb.tip=nb.tip, nb.node=nb.node, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, np=np, f=f, ode.eps=ode.eps, fossil.taxa=fossil.taxa, fix.type=fix.type)
        solution <- numeric(length(pars))
        solution[] <- c(exp(out$solution), 0)[pars]
        
        loglik = -out$objective
    }
    
    names(solution) <- c("turnover00A","turnover01A","turnover10A","turnover11A","eps00A","eps01A","eps10A","eps11A","q00A_01A","q00A_10A","q00A_11A","q01A_00A","q01A_10A","q01A_11A","q10A_00A","q10A_01A","q10A_11A","q11A_00A","q11A_01A","q11A_10A","q00A_00B","q00A_00C","q00A_00D","q00A_00E","q00A_00F","q00A_00G","q00A_00H","q01A_01B","q01A_01C","q01A_01D","q01A_01E","q01A_01F","q01A_01G","q01A_01H","q10A_10B","q10A_10C","q10A_10D","q10A_10E","q10A_10F","q10A_10G","q10A_10H","q11A_11B","q11A_11C","q11A_11D","q11A_11E","q11A_11F","q11A_11G","q11A_11H","turnover00B","turnover01B","turnover10B","turnover11B","eps00B","eps01B","eps10B","eps11B","q00B_01B","q00B_10B","q00B_11B","q01B_00B","q01B_10B","q01B_11B","q10B_00B","q10B_01B","q10B_11B","q11B_00B","q11B_01B","q11B_10B","q00B_00A","q00B_00C","q00B_00D","q00B_00E","q00B_00F","q00B_00G","q00B_00H","q01B_01A","q01B_01C","q01B_01D","q01B_01E","q01B_01F","q01B_01G","q01B_01H","q10B_10A","q10B_10C","q10B_10D","q10B_10E","q10B_10F","q10B_10G","q10B_10H","q11B_11A","q11B_11C","q11B_11D","q11B_11E","q11B_11F","q11B_11G","q11B_11H","turnover00C","turnover01C","turnover10C","turnover11C","eps00C","eps01C","eps10C","eps11C","q00C_01C","q00C_10C","q00C_11C","q01C_00C","q01C_10C","q01C_11C","q10C_00C","q10C_01C","q10C_11C","q11C_00C","q11C_01C","q11C_10C","q00C_00A","q00C_00B","q00C_00D","q00C_00E","q00C_00F","q00C_00G","q00C_00H","q01C_01A","q01C_01B","q01C_01D","q01C_01E","q01C_01F","q01C_01G","q01C_01H","q10C_10A","q10C_10B","q10C_10D","q10C_10E","q10C_10F","q10C_10G","q10C_10H","q11C_11A","q11C_11B","q11C_11D","q11C_11E","q11C_11F","q11C_11G","q11C_11H","turnover00D","turnover01D","turnover10D","turnover11D","eps00D","eps01D","eps10D","eps11D","q00D_01D","q00D_10D","q00D_11D","q01D_00D","q01D_10D","q01D_11D","q10D_00D","q10D_01D","q10D_11D","q11D_00D","q11D_01D","q11D_10D","q00D_00A","q00D_00B","q00D_00C","q00D_00E","q00D_00F","q00D_00G","q00D_00H","q01D_01A","q01D_01B","q01D_01C","q01D_01E","q01D_01F","q01D_01G","q01D_01H","q10D_10A","q10D_10B","q10D_10C","q10D_10E","q10D_10F","q10D_10G","q10D_10H","q11D_11A","q11D_11B","q11D_11C","q11D_11E","q11D_11F","q11D_11G","q11D_11H","turnover00E","turnover01E","turnover10E","turnover11E","eps00E","eps01E","eps10E","eps11E","q00E_01E","q00E_10E","q00E_11E","q01E_00E","q01E_10E","q01E_11E","q10E_00E","q10E_01E","q10E_11E","q11E_00E","q11E_01E","q11E_10E","q00E_00A","q00E_00B","q00E_00C","q00E_00D","q00E_00F","q00E_00G","q00E_00H","q01E_01A","q01E_01B","q01E_01C","q01E_01D","q01E_01F","q01E_01G","q01E_01H","q10E_10A","q10E_10B","q10E_10C","q10E_10D","q10E_10F","q10E_10G","q10E_10H","q11E_11A","q11E_11B","q11E_11C","q11E_11D","q11E_11F","q11E_11G","q11E_11H","turnover00F","turnover01F","turnover10F","turnover11F","eps00F","eps01F","eps10F","eps11F","q00F_01F","q00F_10F","q00F_11F","q01F_00F","q01F_10F","q01F_11F","q10F_00F","q10F_01F","q10F_11F","q11F_00F","q11F_01F","q11F_10F","q00F_00A","q00F_00B","q00F_00C","q00F_00D","q00F_00E","q00F_00G","q00F_00H","q01F_01A","q01F_01B","q01F_01C","q01F_01D","q01F_01E","q01F_01G","q01F_01H","q10F_10A","q10F_10B","q10F_10C","q10F_10D","q10F_10E","q10F_10G","q10F_10H","q11F_11A","q11F_11B","q11F_11C","q11F_11D","q11F_11E","q11F_11G","q11F_11H","turnover00G","turnover01G","turnover10G","turnover11G","eps00G","eps01G","eps10G","eps11G","q00G_01G","q00G_10G","q00G_11G","q01G_00G","q01G_10G","q01G_11G","q10G_00G","q10G_01G","q10G_11G","q11G_00G","q11G_01G","q11G_10G","q00G_00A","q00G_00B","q00G_00C","q00G_00D","q00G_00E","q00G_00F","q00G_00H","q01G_01A","q01G_01B","q01G_01C","q01G_01D","q01G_01E","q01G_01F","q01G_01H","q10G_10A","q10G_10B","q10G_10C","q10G_10D","q10G_10E","q10G_10F","q10G_10H","q11G_11A","q11G_11B","q11G_11C","q11G_11D","q11G_11E","q11G_11F","q11G_11H","turnover00H","turnover01H","turnover10H","turnover11H","eps00H","eps01H","eps10H","eps11H","q00H_01H","q00H_10H","q00H_11H","q01H_00H","q01H_10H","q01H_11H","q10H_00H","q10H_01H","q10H_11H","q11H_00H","q11H_01H","q11H_10H","q00H_00A","q00H_00B","q00H_00C","q00H_00D","q00H_00E","q00H_00F","q00H_00G","q01H_01A","q01H_01B","q01H_01C","q01H_01D","q01H_01E","q01H_01F","q01H_01G","q10H_10A","q10H_10B","q10H_10C","q10H_10D","q10H_10E","q10H_10F","q10H_10G","q11H_11A","q11H_11B","q11H_11C","q11H_11D","q11H_11E","q11H_11F","q11H_11G","psi")
    
    cat("Finished. Summarizing results...", "\n")
    
    obj = list(loglik = loglik, AIC = -2*loglik+2*np, AICc = -2*loglik+(2*np*(Ntip(phy)/(Ntip(phy)-np-1))), solution=solution, index.par=pars, f=f, hidden.states=hidden.states, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, phy=phy, data=data, trans.matrix=trans.rate, max.tol=max.tol, starting.vals=ip, upper.bounds=upper, lower.bounds=lower, ode.eps=ode.eps, includes.fossils=includes.fossils, fix.type=fix.type, psi.type=psi.type, sann)
    class(obj) <- append(class(obj), "muhisse.fit")
    return(obj)
}

######################################################################################################################################
######################################################################################################################################
### The function used to optimize parameters:
######################################################################################################################################
######################################################################################################################################

#Function used for optimizing parameters:
DevOptimizeMuHiSSE <- function(p, pars, dat.tab, gen, hidden.states, nb.tip=nb.tip, nb.node=nb.node, condition.on.survival, root.type, root.p, np, f, ode.eps, fossil.taxa, fix.type) {
    #Generates the final vector with the appropriate parameter estimates in the right place:
    p.new <- exp(p)
    ## print(p.new)
    model.vec <- numeric(length(pars))
    model.vec[] <- c(p.new, 0)[pars]
    cache <- ParametersToPassMuHiSSE(model.vec=model.vec, hidden.states=hidden.states, nb.tip=nb.tip, nb.node=nb.node, bad.likelihood=exp(-300), f=f, ode.eps=ode.eps)
    if(!is.null(fix.type)){
        logl <- DownPassMuHisse(dat.tab=dat.tab, cache=cache, gen=gen, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, node=fix.type[,1], state=NULL, fossil.taxa=fossil.taxa, fix.type=fix.type[,2])
    }else{
        logl <- DownPassMuHisse(dat.tab=dat.tab, cache=cache, gen=gen, condition.on.survival=condition.on.survival, root.type=root.type, root.p=root.p, node=NULL, fossil.taxa=fossil.taxa, fix.type=NULL)
    }
    return(-logl)
}



######################################################################################################################################
######################################################################################################################################
### The various utility functions used
######################################################################################################################################
######################################################################################################################################

OrganizeData <- function(data, phy, f, hidden.states, includes.fossils=FALSE){
    ### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
    DesNode = NULL
    
    nb.tip <- length(phy$tip.label)
    nb.node <- phy$Nnode
    
    if(hidden.states == FALSE){
        states = matrix(0,Ntip(phy),4)
        x <- data[,1]
        y <- data[,2]
        for(i in 1:Ntip(phy)){
            if(x[i]==0 & y[i]==0){states[i,1]=1}
            if(x[i]==0 & y[i]==1){states[i,2]=1}
            if(x[i]==1 & y[i]==0){states[i,3]=1}
            if(x[i]==1 & y[i]==1){states[i,4]=1}
            if(x[i]==2 & y[i]==0){states[i,c(1,3)]=1}
            if(x[i]==2 & y[i]==1){states[i,c(2,4)]=1}
            if(x[i]==0 & y[i]==2){states[i,c(1,2)]=1}
            if(x[i]==1 & y[i]==2){states[i,c(3,4)]=1}
            if(x[i]==2 & y[i]==2){states[i,1:4]=1}
        }
        compD <- matrix(0, nrow=nb.tip, ncol=4)
        compE <- matrix(0, nrow=nb.tip, ncol=4)
    }
    
    if(hidden.states == TRUE){
        states = matrix(0,Ntip(phy),32)
        x <- data[,1]
        y <- data[,2]
        for(i in 1:Ntip(phy)){
            if(x[i]==0 & y[i]==0){states[i,c(1,5,9,13,17,21,25,29)]=1}
            if(x[i]==0 & y[i]==1){states[i,c(2,6,10,14,18,22,26,30)]=1}
            if(x[i]==1 & y[i]==0){states[i,c(3,7,11,15,19,23,27,31)]=1}
            if(x[i]==1 & y[i]==1){states[i,c(4,8,12,16,20,24,28,32)]=1}
            if(x[i]==2 & y[i]==0){states[i,c(1,5,9,13,17,21,25,29, 3,7,11,15,19,23,27,31)]=1}
            if(x[i]==2 & y[i]==1){states[i,c(2,6,10,14,18,22,26,30, 4,8,12,16,20,24,28,32)]=1}
            if(x[i]==0 & y[i]==2){states[i,c(1,5,9,13,17,21,25,29, 2,6,10,14,18,22,26,30)]=1}
            if(x[i]==1 & y[i]==2){states[i,c(3,7,11,15,19,23,27,31, 4,8,12,16,20,24,28,32)]=1}
            if(x[i]==2 & y[i]==2){states[i,1:16]=1}
        }
        compD <- matrix(0, nrow=nb.tip, ncol=32)
        compE <- matrix(0, nrow=nb.tip, ncol=32)
    }
    
    if(hidden.states == "TEST1"){
        states = matrix(0,Ntip(phy),32)
        for(i in 1:Ntip(phy)){
            if(data[i]==1){states[i,c(1,5,9,13,17,21,25,29)]=1}
            if(data[i]==2){states[i,c(2,6,10,14,18,22,26,30)]=1}
            if(data[i]==3){states[i,c(3,7,11,15,19,23,27,31)]=1}
            if(data[i]==4){states[i,c(4,8,12,16,20,24,28,32)]=1}
        }
        compD <- matrix(0, nrow=nb.tip, ncol=32)
        compE <- matrix(0, nrow=nb.tip, ncol=32)
    }
    
    if(hidden.states == "TEST2"){
        states = matrix(0,Ntip(phy),4)
        for(i in 1:Ntip(phy)){
            if(data[i]==1){states[i,1]=1}
            if(data[i]==2){states[i,2]=1}
            if(data[i]==3){states[i,3]=1}
            if(data[i]==4){states[i,4]=1}
        }
        compD <- matrix(0, nrow=nb.tip, ncol=4)
        compE <- matrix(0, nrow=nb.tip, ncol=4)
        
    }
    
    #Initializes the tip sampling and sets internal nodes to be zero:
    ncols = dim(compD)[2]
    if(length(f) == 4){
        for(i in 1:(nb.tip)){
            compD[i,] <- f * states[i,]
            compE[i,] <- rep((1-f), ncols/4)
        }
    }else{
        for(i in 1:(nb.tip)){
            compD[i,] <- f[i] * states[i,]
            compE[i,] <- rep((1-f[i]), ncols/4)
        }
    }
    
    table.info <- GetTreeTable(phy, root.age=NULL)
    if(includes.fossils == TRUE){
        k.sample.tip.no <- grep("Ksamp*", x=phy$tip.label)
        branch.type <- rep(0, dim(table.info)[1])
        for(row.index in 1:dim(table.info)[1]){
            if(table.info[row.index,5]<=nb.tip){
                if(table.info[row.index,2] > .Machine$double.eps^.50){
                    if(includes.fossils == TRUE){
                        branch.type[row.index] <- 1
                    }
                }
                if(any(phy$edge[row.index,2]==k.sample.tip.no)){
                    branch.type[row.index] <- 2
                }
            }
        }
    }else{
        branch.type <- rep(0, dim(table.info)[1])
    }

    tmp.df <- cbind(table.info, 0, matrix(0, nrow(table.info), ncol(compD)), matrix(0, nrow(table.info), ncol(compE)), branch.type)
    colnames(tmp.df) <- c("RootwardAge", "TipwardAge", "BranchLength", "FocalNode", "DesNode", "comp", paste("compD", 1:ncol(compD), sep="_"), paste("compE", 1:ncol(compE), sep="_"), "branch.type")
    dat.tab <- as.data.table(tmp.df)
    setkey(dat.tab, DesNode)
    cols <- names(dat.tab)
    if(hidden.states == TRUE | hidden.states == "TEST1"){
        for (j in 1:(dim(compD)[2])){
            #dat.tab[data.table(c(1:nb.tip)), paste("compD", j, sep="_") := compD[,j]]
            set(dat.tab, 1:nb.tip, cols[6+j], compD[,j])
            #dat.tab[data.table(c(1:nb.tip)), paste("compE", j, sep="_") := compE[,j]]
            set(dat.tab, 1:nb.tip, cols[38+j], compE[,j])
        }
    }else{
        for (j in 1:(dim(compD)[2])){
            #dat.tab[data.table(c(1:nb.tip)), paste("compD", j, sep="_") := compD[,j]]
            set(dat.tab, 1:nb.tip, cols[6+j], compD[,j])
            #dat.tab[data.table(c(1:nb.tip)), paste("compE", j, sep="_") := compE[,j]]
            set(dat.tab, 1:nb.tip, cols[10+j], compE[,j])
        }
    }
    return(dat.tab)
}


SingleChildProb <- function(cache, pars, compD, compE, start.time, end.time, branch.type){
    
    if(any(!is.finite(c(compD, compE)))) { # something went awry at a previous step. Bail!
        prob.subtree.cal <- rep(0, ifelse(cache$hidden.states == TRUE, 64, 8))
        prob.subtree.cal[(1+length(prob.subtree.cal)/2):length(prob.subtree.cal)] <- cache$bad.likelihood
        return(prob.subtree.cal)
    }
    
    if(cache$hidden.states == TRUE){
        yini <- c(E00A = compE[1], E01A = compE[2], E10A = compE[3], E11A = compE[4], E00B = compE[5], E01B = compE[6], E10B = compE[7], E11B = compE[8], E00C = compE[9], E01C = compE[10], E10C = compE[11], E11C = compE[12], E00D = compE[13], E01D = compE[14], E10D = compE[15], E11D = compE[16], E00E = compE[17], E01E = compE[18], E10E = compE[19], E11E = compE[20], E00F = compE[21], E01F = compE[22], E10F = compE[23], E11F = compE[24], E00G = compE[25], E01G = compE[26], E10G = compE[27], E11G = compE[28], E00H = compE[29], E01H = compE[30], E10H = compE[31], E11H = compE[32], D00A = compD[1], D01A = compD[2], D10A = compD[3], D11A = compD[4], D00B = compD[5], D01B = compD[6], D10B = compD[7], D11B = compD[8], D00C = compD[9], D01C = compD[10], D10C = compD[11], D11C = compD[12], D00D = compD[13], D01D = compD[14], D10D = compD[15], D11D = compD[16], D00E = compD[17], D01E = compD[18], D10E = compD[19], D11E = compD[20], D00F = compD[21], D01F = compD[22], D10F = compD[23], D11F = compD[24], D00G = compD[25], D01G = compD[26], D10G = compD[27], D11G = compD[28], D00H = compD[29], D01H = compD[30], D10H = compD[31], D11H = compD[32])
        if(branch.type == 2){
            times <- c(0, end.time)
        }else{
            times=c(start.time, end.time)
        }
        runSilent <- function() {
            options(warn = -1)
            on.exit(options(warn = 0))
            capture.output(res <- lsoda(yini, times, func = "muhisse_derivs", pars, initfunc="initmod_muhisse", dllname = "hisse", rtol=1e-8, atol=1e-8))
            res
        }
        #prob.subtree.cal.full <- lsoda(yini, times, func = "muhisse_derivs", pars, initfunc="initmod_muhisse", dll = "muhisse-ext-derivs", rtol=1e-8, atol=1e-8, hini=10)
        #prob.subtree.cal.full <- lsoda(yini, times, func = "muhisse_derivs", pars, initfunc="initmod_muhisse", dllname = "hisse", rtol=1e-8, atol=1e-8)
        prob.subtree.cal.full <- runSilent()
    }else{
        yini <- c(E00=compE[1], E01=compE[2], E10=compE[3], E11=compE[4], D00=compD[1], D01=compD[2], D10=compD[3], D11=compD[4])
        if(branch.type == 2){
            times <- c(0, end.time)
        }else{
            times=c(start.time, end.time)
        }
        runSilent <- function() {
            options(warn = -1)
            on.exit(options(warn = 0))
            capture.output(res <- lsoda(yini, times, func = "musse_derivs", pars, initfunc="initmod_musse", dllname = "hisse", rtol=1e-8, atol=1e-8))
            res
        }
        #prob.subtree.cal.full <- lsoda(yini, times, func = "musse_derivs", pars, initfunc="initmod_musse", dll = "canonical-musse-ext-derivs", rtol=1e-8, atol=1e-8)
        #prob.subtree.cal.full <- lsoda(yini, times, func = "musse_derivs", pars, initfunc="initmod_musse", dllname = "hisse", rtol=1e-8, atol=1e-8)
        prob.subtree.cal.full <- runSilent()
    }
    
    ######## THIS CHECKS TO ENSURE THAT THE INTEGRATION WAS SUCCESSFUL ###########
    if(attributes(prob.subtree.cal.full)$istate[1] < 0){
        prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
        if(cache$hidden.states == TRUE){
            prob.subtree.cal[33:64] <- cache$bad.likelihood
            return(prob.subtree.cal)
        }else{
            prob.subtree.cal[5:8] <- cache$bad.likelihood
            return(prob.subtree.cal)
        }
    }else{
        if(branch.type == 2){
            prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
            if(cache$hidden.states == TRUE){
                prob.subtree.cal[33:64] <- compD
            }else{
                prob.subtree.cal[5:8] <- compD
            }
        }else{
            if(branch.type == 1){
                if(cache$hidden.states == TRUE){
                    prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
                    prob.subtree.cal[33:64] <- prob.subtree.cal[1:32] * compD
                }else{
                    prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
                    prob.subtree.cal[5:8] <- prob.subtree.cal[1:4] * compD
                }
            }else{
                prob.subtree.cal <- prob.subtree.cal.full[-1,-1]
            }
        }
    }
    ##############################################################################
    
    if(cache$hidden.states == TRUE){
        if(any(is.nan(prob.subtree.cal[33:64]))){
            prob.subtree.cal[33:64] <- cache$bad.likelihood
            return(prob.subtree.cal)
        }
        #This is default and cannot change, but if we get a negative probability, discard the results:
        if(any(prob.subtree.cal[33:64] < 0)){
            prob.subtree.cal[33:64] <- cache$bad.likelihood
            return(prob.subtree.cal)
        }
        if(sum(prob.subtree.cal[33:64]) < cache$ode.eps){
            prob.subtree.cal[33:64] <- cache$bad.likelihood
            return(prob.subtree.cal)
        }
    }else{
        if(any(is.nan(prob.subtree.cal[5:8]))){
            prob.subtree.cal[5:8] <- cache$bad.likelihood
            return(prob.subtree.cal)
        }
        #This is default and cannot change, but if we get a negative probability, discard the results:
        if(any(prob.subtree.cal[5:8] < 0)){
            prob.subtree.cal[5:8] <- cache$bad.likelihood
            return(prob.subtree.cal)
        }
        if(sum(prob.subtree.cal[5:8]) < cache$ode.eps){
            prob.subtree.cal[5:8] <- cache$bad.likelihood
            return(prob.subtree.cal)
        }
    }
    return(prob.subtree.cal)
}


FocalNodeProb <- function(cache, pars, lambdas, dat.tab, generations){
    ### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
    DesNode = NULL
    FocalNode = NULL
    . = NULL
    
    gens <- data.table(c(generations))
    #gens <- dat.tab[.(generations), which=TRUE]
    setkey(dat.tab, FocalNode)
    CurrentGenData <- dat.tab[gens]
    if(cache$hidden.states == TRUE){
        tmp <- t(apply(CurrentGenData, 1, function(z) SingleChildProb(cache, pars, z[7:38], z[39:70], z[2], z[1], z[71])))
        v.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2), 33:64] * tmp[seq(2,nrow(tmp),2), 33:64], length(unique(CurrentGenData$FocalNode)), 32)
        v.mat <- v.mat * matrix(lambdas, length(unique(CurrentGenData$FocalNode)), 32, byrow=TRUE)
        phi.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),1:32], length(unique(CurrentGenData$FocalNode)), 32)
        if(!is.null(cache$node)){
            if(any(cache$node %in% generations)){
                for(fix.index in 1:length(cache$node)){
                    if(cache$fix.type[fix.index] == "event"){
                        #basically we are using the node to fix the state along a branch, but we do not want to assume a true speciation event occurred here.
                        lambdas.check <- lambdas
                        lambdas.check[which(lambdas==0)] <- 1
                        #The initial condition for a k.sample is D(t)*psi
                        v.mat[which(generations == cache$node[fix.index]),] <- (v.mat[which(generations == cache$node[fix.index]),] / lambdas.check) * cache$psi
                    }else{
                        #Fixes the state at the node nothing needs to be done other than fix the node
                        fixer <- numeric(32)
                        fixer[cache$state[fix.index]] = 1
                        v.mat[which(generations == cache$node[fix.index]),] <- v.mat[which(generations == cache$node[fix.index]),] * fixer
                    }
                }
            }
        }
        tmp.comp <- rowSums(v.mat)
        tmp.probs <- v.mat / tmp.comp
        #tmp.probs <- v.mat
        setkey(dat.tab, DesNode)
        #gens <- data.table(c(generations))
        rows <- dat.tab[.(generations), which=TRUE]
        cols <- names(dat.tab)
        for (j in 1:(dim(tmp.probs)[2])){
            #dat.tab[gens, cols[6+j] := tmp.probs[,j]]
            set(dat.tab, rows, cols[6+j], tmp.probs[,j])
            #dat.tab[gens, cols[38+j] := phi.mat[,j]]
            set(dat.tab, rows, cols[38+j], phi.mat[,j])
        }
        dat.tab[gens, "comp" := tmp.comp]
    }else{
        tmp <- t(apply(CurrentGenData, 1, function(z) SingleChildProb(cache, pars, z[7:10], z[11:14], z[2], z[1], z[15])))
        v.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),5:8] * tmp[seq(2,nrow(tmp),2),5:8], length(unique(CurrentGenData$FocalNode)), 4)
        v.mat <- v.mat * matrix(lambdas, length(unique(CurrentGenData$FocalNode)), 4, byrow=TRUE)
        phi.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),1:4], length(unique(CurrentGenData$FocalNode)), 4)
        if(!is.null(cache$node)){
            if(any(cache$node %in% generations)){
                for(fix.index in 1:length(cache$node)){
                    if(cache$fix.type[fix.index] == "event"){
                        #basically we are using the node to fix the state along a branch, but we do not want to assume a true speciation event occurred here.
                        lambdas.check <- lambdas
                        lambdas.check[which(lambdas==0)] <- 1
                        #The initial condition for a k.sample is D(t)*psi
                        v.mat[which(generations == cache$node[fix.index]),] <- (v.mat[which(generations == cache$node[fix.index]),] / lambdas.check) * cache$psi
                    }else{
                        fixer <- numeric(4)
                        fixer[cache$state[fix.index]] = 1
                        #Fixes the state at the node
                        v.mat[which(generations == cache$node[fix.index]),] <- v.mat[which(generations == cache$node[fix.index]),] * fixer
                    }
                }
            }
        }
        tmp.comp <- rowSums(v.mat)
        tmp.probs <- v.mat / tmp.comp
        setkey(dat.tab, DesNode)
        #gens <- data.table(c(generations))
        rows <- dat.tab[.(generations), which=TRUE]
        cols <- names(dat.tab)
        for (j in 1:(dim(tmp.probs)[2])){
            #dat.tab[gens, cols[6+j] := tmp.probs[,j]]
            set(dat.tab, rows, cols[6+j], tmp.probs[,j])
            #dat.tab[gens, cols[10+j] := phi.mat[,j]]
            set(dat.tab, rows, cols[10+j], phi.mat[,j])
        }
        dat.tab[gens, "comp" := tmp.comp]
    }
    return(dat.tab)
}


#Have to calculate root prob separately because it is not a descendant in our table. Could add it, but I worry about the NA that is required.
GetRootProb <- function(cache, pars, lambdas, dat.tab, generations){
    ### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
    FocalNode = NULL
    . = NULL
    
    gens <- data.table(c(generations))
    setkey(dat.tab, FocalNode)
    CurrentGenData <- dat.tab[gens]
    if(cache$hidden.states == TRUE){
        tmp <- t(apply(CurrentGenData, 1, function(z) SingleChildProb(cache, pars, z[7:38], z[39:70], z[2], z[1], z[71])))
        v.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),33:64] * tmp[seq(2,nrow(tmp),2),33:64], length(unique(CurrentGenData$FocalNode)), 32)
        v.mat <- v.mat * matrix(lambdas, length(unique(CurrentGenData$FocalNode)), 32, byrow=TRUE)
        phi.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),1:32], length(unique(CurrentGenData$FocalNode)), 32)
        if(!is.null(cache$node)){
            if(any(cache$node %in% generations)){
                for(fix.index in 1:length(cache$node)){
                    for(fix.index in 1:length(cache$node)){
                        if(cache$fix.type[fix.index] == "event"){
                            #basically we are using the node to fix the state along a branch, but we do not want to assume a true speciation event occurred here.
                            lambdas.check <- lambdas
                            lambdas.check[which(lambdas==0)] <- 1
                            #The initial condition for a k.sample is D(t)*psi
                            v.mat[which(generations == cache$node[fix.index]),] <- (v.mat[which(generations == cache$node[fix.index]),] / lambdas.check) * cache$psi
                        }else{
                            #Fixes the state at the node nothing needs to be done other than fix the node
                            fixer = numeric(32)
                            fixer[cache$state[fix.index]] = 1
                            v.mat[which(generations == cache$node[fix.index]),] <- v.mat[which(generations == cache$node[fix.index]),] * fixer
                        }
                    }
                }
            }
        }
    }else{
        tmp <- t(apply(CurrentGenData, 1, function(z) SingleChildProb(cache, pars, z[7:10], z[11:14], z[2], z[1], z[15])))
        v.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),5:8] * tmp[seq(2,nrow(tmp),2),5:8], length(unique(CurrentGenData$FocalNode)), 4)
        v.mat <- v.mat * matrix(lambdas, length(unique(CurrentGenData$FocalNode)), 4, byrow=TRUE)
        phi.mat <- matrix(tmp[seq(1,nrow(tmp)-1,2),1:4], length(unique(CurrentGenData$FocalNode)), 4)
        if(!is.null(cache$node)){
            if(any(cache$node %in% generations)){
                for(fix.index in 1:length(cache$node)){
                    if(cache$fix.type[fix.index] == "event"){
                        lambdas.check <- lambdas
                        lambdas.check[which(lambdas==0)] <- 1
                        #The initial condition for a k.sample is D(t)*psi
                        v.mat[which(generations == cache$node[fix.index]),] <- (v.mat[which(generations == cache$node[fix.index]),] / lambdas.check) * cache$psi
                    }else{
                        fixer = numeric(4)
                        fixer[cache$state[fix.index]] = 1
                        #Fixes the state at the node
                        v.mat[which(generations == cache$node[fix.index]),] <- v.mat[which(generations == cache$node[fix.index]),] * fixer
                    }
                }
            }
        }
    }
    
    tmp.comp <- rowSums(v.mat)
    tmp.probs <- v.mat / tmp.comp
    #tmp.probs <- v.mat
    
    return(cbind(tmp.comp, phi.mat, tmp.probs))
}


#Calculates the initial conditions for fossil taxa in the tree.
GetFossilInitialsMuHiSSE <- function(cache, pars, lambdas, dat.tab, fossil.taxa){
    ### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
    DesNode = NULL
    . = NULL
    if(cache$hidden.states == TRUE){
        fossils <- data.table(c(fossil.taxa))
        setkey(dat.tab, DesNode)
        CurrentGenData <- dat.tab[fossils]
        tmp <- t(apply(CurrentGenData, 1, function(z) SingleChildProb(cache, pars, z[7:38], z[39:70],  0, z[2], z[71])))
        tmp.probs <- matrix(tmp[,33:64], length(fossil.taxa), 32) * cache$psi
        phi.mat <- matrix(tmp[,1:32], length(fossil.taxa), 32)
        
        setkey(dat.tab, DesNode)
        rows <- dat.tab[.(fossils), which=TRUE]
        cols <- names(dat.tab)
        for (j in 1:(dim(tmp.probs)[2])){
            #dat.tab[data.table(c(generations)), paste("compD", j, sep="_") := tmp.probs[,j]]
            set(dat.tab, rows, cols[6+j], tmp.probs[,j])
            #dat.tab[data.table(c(generations)), paste("compE", j, sep="_") := phi.mat[,j]]
            set(dat.tab, rows, cols[38+j], phi.mat[,j])
        }
    }else{
        fossils <- data.table(c(fossil.taxa))
        setkey(dat.tab, DesNode)
        CurrentGenData <- dat.tab[fossils]
        tmp <- t(apply(CurrentGenData, 1, function(z) SingleChildProb(cache, pars, z[7:10], z[11:14], 0, z[2], z[15])))
        tmp.probs <- matrix(tmp[,5:8], length(fossil.taxa), 4) * cache$psi
        phi.mat <- matrix(tmp[,1:4], length(fossil.taxa), 4)
        
        setkey(dat.tab, DesNode)
        rows <- dat.tab[.(fossils), which=TRUE]
        cols <- names(dat.tab)
        for (j in 1:(dim(tmp.probs)[2])){
            #dat.tab[data.table(c(generations)), paste("compD", j, sep="_") := tmp.probs[,j]]
            set(dat.tab, rows, cols[6+j], tmp.probs[,j])
            #dat.tab[gens, cols[10+j] := phi.mat[,j]]
            set(dat.tab, rows, cols[10+j], phi.mat[,j])
        }
    }
    dat.tab[fossils,"branch.type" := 0]
    return(dat.tab)
}



######################################################################################################################################
######################################################################################################################################
### The MuHiSSE type down pass that carries out the integration and returns the likelihood:
######################################################################################################################################
######################################################################################################################################

DownPassMuHisse <- function(dat.tab, gen, cache, condition.on.survival, root.type, root.p, get.phi=FALSE, node=NULL, state=NULL, fossil.taxa=NULL, fix.type=NULL) {
    
    ### Ughy McUgherson. This is a must in order to pass CRAN checks: http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
    DesNode = NULL
    compE = NULL
    . = NULL
    
    ## Moved this here instead of above because it actually significantly improved the speed.
    if(cache$hidden.states == FALSE){
        pars <- c(cache$lambda00A, cache$lambda01A, cache$lambda10A, cache$lambda11A, cache$mu00A, cache$mu01A, cache$mu10A, cache$mu11A, cache$q00A_01A, cache$q00A_10A, cache$q00A_11A, cache$q01A_00A, cache$q01A_10A, cache$q01A_11A, cache$q10A_00A, cache$q10A_01A, cache$q10A_11A, cache$q11A_00A, cache$q11A_01A, cache$q11A_10A, cache$psi)
        lambda <- c(cache$lambda00A, cache$lambda01A, cache$lambda10A, cache$lambda11A)
    }else{
        pars <- c(cache$lambda00A, cache$lambda01A, cache$lambda10A, cache$lambda11A, cache$mu00A, cache$mu01A, cache$mu10A, cache$mu11A, cache$q00A_01A, cache$q00A_10A, cache$q00A_11A, cache$q01A_00A, cache$q01A_10A, cache$q01A_11A, cache$q10A_00A, cache$q10A_01A, cache$q10A_11A, cache$q11A_00A, cache$q11A_01A, cache$q11A_10A, cache$q00A_00B, cache$q00A_00C, cache$q00A_00D, cache$q00A_00E, cache$q00A_00F, cache$q00A_00G, cache$q00A_00H, cache$q01A_01B, cache$q01A_01C, cache$q01A_01D, cache$q01A_01E, cache$q01A_01F, cache$q01A_01G, cache$q01A_01H, cache$q10A_10B, cache$q10A_10C, cache$q10A_10D, cache$q10A_10E, cache$q10A_10F, cache$q10A_10G, cache$q10A_10H, cache$q11A_11B, cache$q11A_11C, cache$q11A_11D, cache$q11A_11E, cache$q11A_11F, cache$q11A_11G, cache$q11A_11H, cache$lambda00B, cache$lambda01B, cache$lambda10B, cache$lambda11B, cache$mu00B, cache$mu01B, cache$mu10B, cache$mu11B, cache$q00B_01B, cache$q00B_10B, cache$q00B_11B, cache$q01B_00B, cache$q01B_10B, cache$q01B_11B, cache$q10B_00B, cache$q10B_01B, cache$q10B_11B, cache$q11B_00B, cache$q11B_01B, cache$q11B_10B, cache$q00B_00A, cache$q00B_00C, cache$q00B_00D, cache$q00B_00E, cache$q00B_00F, cache$q00B_00G, cache$q00B_00H, cache$q01B_01A, cache$q01B_01C, cache$q01B_01D, cache$q01B_01E, cache$q01B_01F, cache$q01B_01G, cache$q01B_01H, cache$q10B_10A, cache$q10B_10C, cache$q10B_10D, cache$q10B_10E, cache$q10B_10F, cache$q10B_10G, cache$q10B_10H, cache$q11B_11A, cache$q11B_11C, cache$q11B_11D, cache$q11B_11E, cache$q11B_11F, cache$q11B_11G, cache$q11B_11H, cache$lambda00C, cache$lambda01C, cache$lambda10C, cache$lambda11C, cache$mu00C, cache$mu01C, cache$mu10C, cache$mu11C, cache$q00C_01C, cache$q00C_10C, cache$q00C_11C, cache$q01C_00C, cache$q01C_10C, cache$q01C_11C, cache$q10C_00C, cache$q10C_01C, cache$q10C_11C, cache$q11C_00C, cache$q11C_01C, cache$q11C_10C, cache$q00C_00A, cache$q00C_00B, cache$q00C_00D, cache$q00C_00E, cache$q00C_00F, cache$q00C_00G, cache$q00C_00H, cache$q01C_01A, cache$q01C_01B, cache$q01C_01D, cache$q01C_01E, cache$q01C_01F, cache$q01C_01G, cache$q01C_01H, cache$q10C_10A, cache$q10C_10B, cache$q10C_10D, cache$q10C_10E, cache$q10C_10F, cache$q10C_10G, cache$q10C_10H, cache$q11C_11A, cache$q11C_11B, cache$q11C_11D, cache$q11C_11E, cache$q11C_11F, cache$q11C_11G, cache$q11C_11H, cache$lambda00D, cache$lambda01D, cache$lambda10D, cache$lambda11D, cache$mu00D, cache$mu01D, cache$mu10D, cache$mu11D, cache$q00D_01D, cache$q00D_10D, cache$q00D_11D, cache$q01D_00D, cache$q01D_10D, cache$q01D_11D, cache$q10D_00D, cache$q10D_01D, cache$q10D_11D, cache$q11D_00D, cache$q11D_01D, cache$q11D_10D, cache$q00D_00A, cache$q00D_00B, cache$q00D_00C, cache$q00D_00E, cache$q00D_00F, cache$q00D_00G, cache$q00D_00H, cache$q01D_01A, cache$q01D_01B, cache$q01D_01C, cache$q01D_01E, cache$q01D_01F, cache$q01D_01G, cache$q01D_01H, cache$q10D_10A, cache$q10D_10B, cache$q10D_10C, cache$q10D_10E, cache$q10D_10F, cache$q10D_10G, cache$q10D_10H, cache$q11D_11A, cache$q11D_11B, cache$q11D_11C, cache$q11D_11E, cache$q11D_11F, cache$q11D_11G, cache$q11D_11H, cache$lambda00E, cache$lambda01E, cache$lambda10E, cache$lambda11E, cache$mu00E, cache$mu01E, cache$mu10E, cache$mu11E, cache$q00E_01E, cache$q00E_10E, cache$q00E_11E, cache$q01E_00E, cache$q01E_10E, cache$q01E_11E, cache$q10E_00E, cache$q10E_01E, cache$q10E_11E, cache$q11E_00E, cache$q11E_01E, cache$q11E_10E, cache$q00E_00A, cache$q00E_00B, cache$q00E_00C, cache$q00E_00D, cache$q00E_00F, cache$q00E_00G, cache$q00E_00H, cache$q01E_01A, cache$q01E_01B, cache$q01E_01C, cache$q01E_01D, cache$q01E_01F, cache$q01E_01G, cache$q01E_01H, cache$q10E_10A, cache$q10E_10B, cache$q10E_10C, cache$q10E_10D, cache$q10E_10F, cache$q10E_10G, cache$q10E_10H, cache$q11E_11A, cache$q11E_11B, cache$q11E_11C, cache$q11E_11D, cache$q11E_11F, cache$q11E_11G, cache$q11E_11H, cache$lambda00F, cache$lambda01F, cache$lambda10F, cache$lambda11F, cache$mu00F, cache$mu01F, cache$mu10F, cache$mu11F, cache$q00F_01F, cache$q00F_10F, cache$q00F_11F, cache$q01F_00F, cache$q01F_10F, cache$q01F_11F, cache$q10F_00F, cache$q10F_01F, cache$q10F_11F, cache$q11F_00F, cache$q11F_01F, cache$q11F_10F, cache$q00F_00A, cache$q00F_00B, cache$q00F_00C, cache$q00F_00D, cache$q00F_00E, cache$q00F_00G, cache$q00F_00H, cache$q01F_01A, cache$q01F_01B, cache$q01F_01C, cache$q01F_01D, cache$q01F_01E, cache$q01F_01G, cache$q01F_01H, cache$q10F_10A, cache$q10F_10B, cache$q10F_10C, cache$q10F_10D, cache$q10F_10E, cache$q10F_10G, cache$q10F_10H, cache$q11F_11A, cache$q11F_11B, cache$q11F_11C, cache$q11F_11D, cache$q11F_11E, cache$q11F_11G, cache$q11F_11H, cache$lambda00G, cache$lambda01G, cache$lambda10G, cache$lambda11G, cache$mu00G, cache$mu01G, cache$mu10G, cache$mu11G, cache$q00G_01G, cache$q00G_10G, cache$q00G_11G, cache$q01G_00G, cache$q01G_10G, cache$q01G_11G, cache$q10G_00G, cache$q10G_01G, cache$q10G_11G, cache$q11G_00G, cache$q11G_01G, cache$q11G_10G, cache$q00G_00A, cache$q00G_00B, cache$q00G_00C, cache$q00G_00D, cache$q00G_00E, cache$q00G_00F, cache$q00G_00H, cache$q01G_01A, cache$q01G_01B, cache$q01G_01C, cache$q01G_01D, cache$q01G_01E, cache$q01G_01F, cache$q01G_01H, cache$q10G_10A, cache$q10G_10B, cache$q10G_10C, cache$q10G_10D, cache$q10G_10E, cache$q10G_10F, cache$q10G_10H, cache$q11G_11A, cache$q11G_11B, cache$q11G_11C, cache$q11G_11D, cache$q11G_11E, cache$q11G_11F, cache$q11G_11H, cache$lambda00H, cache$lambda01H, cache$lambda10H, cache$lambda11H, cache$mu00H, cache$mu01H, cache$mu10H, cache$mu11H, cache$q00H_01H, cache$q00H_10H, cache$q00H_11H, cache$q01H_00H, cache$q01H_10H, cache$q01H_11H, cache$q10H_00H, cache$q10H_01H, cache$q10H_11H, cache$q11H_00H, cache$q11H_01H, cache$q11H_10H, cache$q00H_00A, cache$q00H_00B, cache$q00H_00C, cache$q00H_00D, cache$q00H_00E, cache$q00H_00F, cache$q00H_00G, cache$q01H_01A, cache$q01H_01B, cache$q01H_01C, cache$q01H_01D, cache$q01H_01E, cache$q01H_01F, cache$q01H_01G, cache$q10H_10A, cache$q10H_10B, cache$q10H_10C, cache$q10H_10D, cache$q10H_10E, cache$q10H_10F, cache$q10H_10G, cache$q11H_11A, cache$q11H_11B, cache$q11H_11C, cache$q11H_11D, cache$q11H_11E, cache$q11H_11F, cache$q11H_11G, cache$psi)
        lambda <- c(cache$lambda00A,cache$lambda01A,cache$lambda10A,cache$lambda11A,cache$lambda00B,cache$lambda01B,cache$lambda10B,cache$lambda11B,cache$lambda00C,cache$lambda01C,cache$lambda10C,cache$lambda11C,cache$lambda00D,cache$lambda01D,cache$lambda10D,cache$lambda11D,cache$lambda00E,cache$lambda01E,cache$lambda10E,cache$lambda11E,cache$lambda00F,cache$lambda01F,cache$lambda10F,cache$lambda11F,cache$lambda00G,cache$lambda01G,cache$lambda10G,cache$lambda11G,cache$lambda00H,cache$lambda01H,cache$lambda10H,cache$lambda11H)
    }
    
    if(!is.null(fossil.taxa)){
        #Gets initial conditions for fossil taxa:
        dat.tab.copy <- copy(dat.tab)
        dat.tab.copy <- GetFossilInitialsMuHiSSE(cache=cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, fossil.taxa=fossil.taxa)
    }else{
        dat.tab.copy <- copy(dat.tab)
    }
    
    TIPS <- 1:cache$nb.tip
    for(i in 1:length(gen)){
        if(i == length(gen)){
            if(!is.null(node)){
                if(any(node %in% gen[[i]])){
                    cache$node <- node
                    cache$state <- state
                    cache$fix.type <- fix.type
                    res.tmp <- GetRootProb(cache=cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, generations=gen[[i]])
                    cache$node <- NULL
                    cache$state <- NULL
                    cache$fix.type <- NULL
                }else{
                    res.tmp <- GetRootProb(cache=cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, generations=gen[[i]])
                }
            }else{
                res.tmp <- GetRootProb(cache=cache, pars=pars, lambdas=lambda, dat.tab=dat.tab.copy, generations=gen[[i]])
            }
            if(cache$hidden.states == TRUE){
                compD.root <- res.tmp[c(34:65)]
                compE.root <- res.tmp[c(2:33)]
            }else{
                compD.root <- res.tmp[c(6:9)]
                compE.root <- res.tmp[c(2:5)]
            }
            setkey(dat.tab.copy, DesNode)
            comp <- dat.tab.copy[["comp"]]
            comp <- c(comp[-TIPS], res.tmp[1])
        }else{
            if(!is.null(node)){
                if(any(node %in% gen[[i]])){
                    cache$node <- node
                    cache$state <- state
                    cache$fix.type <- fix.type
                    dat.tab.copy <- FocalNodeProb(cache, pars=pars, lambdas=lambda, dat.tab.copy, gen[[i]])
                    cache$node <- NULL
                    cache$state <- NULL
                    cache$fix.type <- NULL
                }else{
                    dat.tab.copy <- FocalNodeProb(cache, pars=pars, lambdas=lambda, dat.tab.copy, gen[[i]])
                }
            }else{
                dat.tab.copy <- FocalNodeProb(cache, pars=pars, lambdas=lambda, dat.tab.copy, gen[[i]])
            }
        }
    }
    
    if(!is.null(fossil.taxa)){
        if(cache$hidden.states == TRUE){
            pars[length(pars)] <- 0
            cache$psi <- 0
            init.d <- rep(cache$f, 8)
            init.e <- rep(1-cache$f, 8)
            phi.mat <- SingleChildProb(cache, pars, init.d, init.e, 0, max(dat.tab$RootwardAge), 0)
            compE.root <- matrix(phi.mat[1:32], 1, 32)
        }else{
            pars[length(pars)] <- 0
            cache$psi <- 0
            init.d <- cache$f
            init.e <- 1-cache$f
            phi.mat <- SingleChildProb(cache, pars, init.d, init.e, 0, max(dat.tab$RootwardAge), 0)
            compE.root <- matrix(phi.mat[1:4], 1, 4)
        }
    }
    
    if (is.na(sum(log(compD.root))) || is.na(log(sum(1-compE.root)))){
        return(log(cache$bad.likelihood)^13)
    }else{
        if(root.type == "madfitz" | root.type == "herr_als"){
            if(is.null(root.p)){
                root.p = compD.root/sum(compD.root)
                root.p[which(is.na(root.p))] = 0
            }
        }
        if(condition.on.survival == TRUE){
            if(cache$hidden.states == FALSE){
                if(root.type == "madfitz"){
                    compD.root <- compD.root / sum(root.p * lambda * (1 - compE.root)^2)
                    #Corrects for possibility that you have 0/0:
                    compD.root[which(is.na(compD.root))] = 0
                    loglik <- log(sum(compD.root * root.p)) + sum(log(comp))
                }else{
                    compD.root <- (compD.root * root.p) / (lambda * (1 - compE.root)^2)
                    #Corrects for possibility that you have 0/0:
                    compD.root[which(is.na(compD.root))] = 0
                    loglik <- log(sum(compD.root)) + sum(log(comp))
                }
            }else{
                if(root.type == "madfitz"){
                    compD.root <- compD.root / sum(root.p * lambda * (1 - compE.root)^2)
                    #Corrects for possibility that you have 0/0:
                    compD.root[which(is.na(compD.root))] = 0
                    loglik <- log(sum(compD.root * root.p)) + sum(log(comp))
                    #Here for testing purposes:
                    #loglik <- log(sum(compD.root * root.p))
                }else{
                    compD.root <- (compD.root * root.p) / (lambda * (1 - compE.root)^2)
                    #Corrects for possibility that you have 0/0:
                    compD.root[which(is.na(compD.root))] = 0
                    loglik <- log(sum(compD.root)) + sum(log(comp))
                }
            }
        }
        if(!is.finite(loglik)){
            return(log(cache$bad.likelihood)^7)
        }
    }
    if(get.phi==TRUE){
        obj = NULL
        obj$compD.root = compD.root/sum(compD.root)
        obj$compE = compE.root
        obj$root.p = root.p
        return(obj)
    }else{
        return(loglik)
    }
}


######################################################################################################################################
######################################################################################################################################
### Cache object for storing parameters that are used throughout MuHiSSE:
######################################################################################################################################
######################################################################################################################################

ParametersToPassMuHiSSE <- function(model.vec, hidden.states, nb.tip, nb.node, bad.likelihood, f, ode.eps){
    #Provides an initial object that contains all the parameters to be passed among functions. This will also be used to pass other things are we move down the tree (see DownPassGeoSSE):
    obj <- NULL
    
    obj$hidden.states <- hidden.states
    obj$nb.tip <- nb.tip
    obj$nb.node <- nb.node
    obj$bad.likelihood <- bad.likelihood
    obj$ode.eps <- ode.eps
    obj$f <- f
    
    ##Hidden State A
    obj$lambda00A = model.vec[1] / (1 + model.vec[5])
    obj$lambda01A = model.vec[2] / (1 + model.vec[6])
    obj$lambda10A = model.vec[3] / (1 + model.vec[7])
    obj$lambda11A = model.vec[4] / (1 + model.vec[8])
    obj$mu00A = (model.vec[5] * model.vec[1]) / (1 + model.vec[5])
    obj$mu01A = (model.vec[6] * model.vec[2]) / (1 + model.vec[6])
    obj$mu10A = (model.vec[7] * model.vec[3]) / (1 + model.vec[7])
    obj$mu11A = (model.vec[8] * model.vec[4]) / (1 + model.vec[8])
    obj$q00A_01A = model.vec[9]
    obj$q00A_10A = model.vec[10]
    obj$q00A_11A = model.vec[11]
    obj$q01A_00A = model.vec[12]
    obj$q01A_10A = model.vec[13]
    obj$q01A_11A = model.vec[14]
    obj$q10A_00A = model.vec[15]
    obj$q10A_01A = model.vec[16]
    obj$q10A_11A = model.vec[17]
    obj$q11A_00A = model.vec[18]
    obj$q11A_01A = model.vec[19]
    obj$q11A_10A = model.vec[20]
    
    obj$q00A_00B = model.vec[21]
    obj$q00A_00C = model.vec[22]
    obj$q00A_00D = model.vec[23]
    obj$q00A_00E = model.vec[24]
    obj$q00A_00F = model.vec[25]
    obj$q00A_00G = model.vec[26]
    obj$q00A_00H = model.vec[27]
    obj$q01A_01B = model.vec[28]
    obj$q01A_01C = model.vec[29]
    obj$q01A_01D = model.vec[30]
    obj$q01A_01E = model.vec[31]
    obj$q01A_01F = model.vec[32]
    obj$q01A_01G = model.vec[33]
    obj$q01A_01H = model.vec[34]
    obj$q10A_10B = model.vec[35]
    obj$q10A_10C = model.vec[36]
    obj$q10A_10D = model.vec[37]
    obj$q10A_10E = model.vec[38]
    obj$q10A_10F = model.vec[39]
    obj$q10A_10G = model.vec[40]
    obj$q10A_10H = model.vec[41]
    obj$q11A_11B = model.vec[42]
    obj$q11A_11C = model.vec[43]
    obj$q11A_11D = model.vec[44]
    obj$q11A_11E = model.vec[45]
    obj$q11A_11F = model.vec[46]
    obj$q11A_11G = model.vec[47]
    obj$q11A_11H = model.vec[48]
    
    ##Hidden State B
    obj$lambda00B = model.vec[49] / (1 + model.vec[53])
    obj$lambda01B = model.vec[50] / (1 + model.vec[54])
    obj$lambda10B = model.vec[51] / (1 + model.vec[55])
    obj$lambda11B = model.vec[52] / (1 + model.vec[56])
    obj$mu00B = (model.vec[53] * model.vec[49]) / (1 + model.vec[53])
    obj$mu01B = (model.vec[54] * model.vec[50]) / (1 + model.vec[54])
    obj$mu10B = (model.vec[55] * model.vec[51]) / (1 + model.vec[55])
    obj$mu11B = (model.vec[56] * model.vec[52]) / (1 + model.vec[56])
    obj$q00B_01B = model.vec[57]
    obj$q00B_10B = model.vec[58]
    obj$q00B_11B = model.vec[59]
    obj$q01B_00B = model.vec[60]
    obj$q01B_10B = model.vec[61]
    obj$q01B_11B = model.vec[62]
    obj$q10B_00B = model.vec[63]
    obj$q10B_01B = model.vec[64]
    obj$q10B_11B = model.vec[65]
    obj$q11B_00B = model.vec[66]
    obj$q11B_01B = model.vec[67]
    obj$q11B_10B = model.vec[68]
    
    obj$q00B_00A = model.vec[69]
    obj$q00B_00C = model.vec[70]
    obj$q00B_00D = model.vec[71]
    obj$q00B_00E = model.vec[72]
    obj$q00B_00F = model.vec[73]
    obj$q00B_00G = model.vec[74]
    obj$q00B_00H = model.vec[75]
    obj$q01B_01A = model.vec[76]
    obj$q01B_01C = model.vec[77]
    obj$q01B_01D = model.vec[78]
    obj$q01B_01E = model.vec[79]
    obj$q01B_01F = model.vec[80]
    obj$q01B_01G = model.vec[81]
    obj$q01B_01H = model.vec[82]
    obj$q10B_10A = model.vec[83]
    obj$q10B_10C = model.vec[84]
    obj$q10B_10D = model.vec[85]
    obj$q10B_10E = model.vec[86]
    obj$q10B_10F = model.vec[87]
    obj$q10B_10G = model.vec[88]
    obj$q10B_10H = model.vec[89]
    obj$q11B_11A = model.vec[90]
    obj$q11B_11C = model.vec[91]
    obj$q11B_11D = model.vec[92]
    obj$q11B_11E = model.vec[93]
    obj$q11B_11F = model.vec[94]
    obj$q11B_11G = model.vec[95]
    obj$q11B_11H = model.vec[96]
    
    ##Hidden State C
    obj$lambda00C = model.vec[97] / (1 + model.vec[101])
    obj$lambda01C = model.vec[98] / (1 + model.vec[102])
    obj$lambda10C = model.vec[99] / (1 + model.vec[103])
    obj$lambda11C = model.vec[100] / (1 + model.vec[104])
    obj$mu00C = (model.vec[101] * model.vec[97]) / (1 + model.vec[101])
    obj$mu01C = (model.vec[102] * model.vec[98]) / (1 + model.vec[102])
    obj$mu10C = (model.vec[103] * model.vec[99]) / (1 + model.vec[103])
    obj$mu11C = (model.vec[104] * model.vec[100]) / (1 + model.vec[104])
    obj$q00C_01C = model.vec[105]
    obj$q00C_10C = model.vec[106]
    obj$q00C_11C = model.vec[107]
    obj$q01C_00C = model.vec[108]
    obj$q01C_10C = model.vec[109]
    obj$q01C_11C = model.vec[110]
    obj$q10C_00C = model.vec[111]
    obj$q10C_01C = model.vec[112]
    obj$q10C_11C = model.vec[113]
    obj$q11C_00C = model.vec[114]
    obj$q11C_01C = model.vec[115]
    obj$q11C_10C = model.vec[116]
    
    obj$q00C_00A = model.vec[117]
    obj$q00C_00B = model.vec[118]
    obj$q00C_00D = model.vec[119]
    obj$q00C_00E = model.vec[120]
    obj$q00C_00F = model.vec[121]
    obj$q00C_00G = model.vec[122]
    obj$q00C_00H = model.vec[123]
    obj$q01C_01A = model.vec[124]
    obj$q01C_01B = model.vec[125]
    obj$q01C_01D = model.vec[126]
    obj$q01C_01E = model.vec[127]
    obj$q01C_01F = model.vec[128]
    obj$q01C_01G = model.vec[129]
    obj$q01C_01H = model.vec[130]
    obj$q10C_10A = model.vec[131]
    obj$q10C_10B = model.vec[132]
    obj$q10C_10D = model.vec[133]
    obj$q10C_10E = model.vec[134]
    obj$q10C_10F = model.vec[135]
    obj$q10C_10G = model.vec[136]
    obj$q10C_10H = model.vec[137]
    obj$q11C_11A = model.vec[138]
    obj$q11C_11B = model.vec[139]
    obj$q11C_11D = model.vec[140]
    obj$q11C_11E = model.vec[141]
    obj$q11C_11F = model.vec[142]
    obj$q11C_11G = model.vec[143]
    obj$q11C_11H = model.vec[144]
    
    ##Hidden State D
    obj$lambda00D = model.vec[145] / (1 + model.vec[149])
    obj$lambda01D = model.vec[146] / (1 + model.vec[150])
    obj$lambda10D = model.vec[147] / (1 + model.vec[151])
    obj$lambda11D = model.vec[148] / (1 + model.vec[152])
    obj$mu00D = (model.vec[149] * model.vec[145]) / (1 + model.vec[149])
    obj$mu01D = (model.vec[150] * model.vec[146]) / (1 + model.vec[150])
    obj$mu10D = (model.vec[151] * model.vec[147]) / (1 + model.vec[151])
    obj$mu11D = (model.vec[152] * model.vec[148]) / (1 + model.vec[152])
    obj$q00D_01D = model.vec[153]
    obj$q00D_10D = model.vec[154]
    obj$q00D_11D = model.vec[155]
    obj$q01D_00D = model.vec[156]
    obj$q01D_10D = model.vec[157]
    obj$q01D_11D = model.vec[158]
    obj$q10D_00D = model.vec[159]
    obj$q10D_01D = model.vec[160]
    obj$q10D_11D = model.vec[161]
    obj$q11D_00D = model.vec[162]
    obj$q11D_01D = model.vec[163]
    obj$q11D_10D = model.vec[164]
    
    obj$q00D_00A = model.vec[165]
    obj$q00D_00B = model.vec[166]
    obj$q00D_00C = model.vec[167]
    obj$q00D_00E = model.vec[168]
    obj$q00D_00F = model.vec[169]
    obj$q00D_00G = model.vec[170]
    obj$q00D_00H = model.vec[171]
    obj$q01D_01A = model.vec[172]
    obj$q01D_01B = model.vec[173]
    obj$q01D_01C = model.vec[174]
    obj$q01D_01E = model.vec[175]
    obj$q01D_01F = model.vec[176]
    obj$q01D_01G = model.vec[177]
    obj$q01D_01H = model.vec[178]
    obj$q10D_10A = model.vec[179]
    obj$q10D_10B = model.vec[180]
    obj$q10D_10C = model.vec[181]
    obj$q10D_10E = model.vec[182]
    obj$q10D_10F = model.vec[183]
    obj$q10D_10G = model.vec[184]
    obj$q10D_10H = model.vec[185]
    obj$q11D_11A = model.vec[186]
    obj$q11D_11B = model.vec[187]
    obj$q11D_11C = model.vec[188]
    obj$q11D_11E = model.vec[189]
    obj$q11D_11F = model.vec[190]
    obj$q11D_11G = model.vec[191]
    obj$q11D_11H = model.vec[192]
    
    ##Hidden State E
    obj$lambda00E = model.vec[193] / (1 + model.vec[197])
    obj$lambda01E = model.vec[194] / (1 + model.vec[198])
    obj$lambda10E = model.vec[195] / (1 + model.vec[199])
    obj$lambda11E = model.vec[196] / (1 + model.vec[200])
    obj$mu00E = (model.vec[197] * model.vec[193]) / (1 + model.vec[197])
    obj$mu01E = (model.vec[198] * model.vec[194]) / (1 + model.vec[198])
    obj$mu10E = (model.vec[199] * model.vec[195]) / (1 + model.vec[199])
    obj$mu11E = (model.vec[200] * model.vec[196]) / (1 + model.vec[200])
    obj$q00E_01E = model.vec[201]
    obj$q00E_10E = model.vec[202]
    obj$q00E_11E = model.vec[203]
    obj$q01E_00E = model.vec[204]
    obj$q01E_10E = model.vec[205]
    obj$q01E_11E = model.vec[206]
    obj$q10E_00E = model.vec[207]
    obj$q10E_01E = model.vec[208]
    obj$q10E_11E = model.vec[209]
    obj$q11E_00E = model.vec[210]
    obj$q11E_01E = model.vec[211]
    obj$q11E_10E = model.vec[212]
    
    obj$q00E_00A = model.vec[213]
    obj$q00E_00B = model.vec[214]
    obj$q00E_00C = model.vec[215]
    obj$q00E_00D = model.vec[216]
    obj$q00E_00F = model.vec[217]
    obj$q00E_00G = model.vec[218]
    obj$q00E_00H = model.vec[219]
    obj$q01E_01A = model.vec[220]
    obj$q01E_01B = model.vec[221]
    obj$q01E_01C = model.vec[222]
    obj$q01E_01D = model.vec[223]
    obj$q01E_01F = model.vec[224]
    obj$q01E_01G = model.vec[225]
    obj$q01E_01H = model.vec[226]
    obj$q10E_10A = model.vec[227]
    obj$q10E_10B = model.vec[228]
    obj$q10E_10C = model.vec[229]
    obj$q10E_10D = model.vec[230]
    obj$q10E_10F = model.vec[231]
    obj$q10E_10G = model.vec[232]
    obj$q10E_10H = model.vec[233]
    obj$q11E_11A = model.vec[234]
    obj$q11E_11B = model.vec[235]
    obj$q11E_11C = model.vec[236]
    obj$q11E_11D = model.vec[237]
    obj$q11E_11F = model.vec[238]
    obj$q11E_11G = model.vec[239]
    obj$q11E_11H = model.vec[240]
    
    ##Hidden State F
    obj$lambda00F = model.vec[241] / (1 + model.vec[245])
    obj$lambda01F = model.vec[242] / (1 + model.vec[246])
    obj$lambda10F = model.vec[243] / (1 + model.vec[247])
    obj$lambda11F = model.vec[244] / (1 + model.vec[248])
    obj$mu00F = (model.vec[245] * model.vec[241]) / (1 + model.vec[245])
    obj$mu01F = (model.vec[246] * model.vec[242]) / (1 + model.vec[246])
    obj$mu10F = (model.vec[247] * model.vec[243]) / (1 + model.vec[247])
    obj$mu11F = (model.vec[248] * model.vec[244]) / (1 + model.vec[248])
    obj$q00F_01F = model.vec[249]
    obj$q00F_10F = model.vec[250]
    obj$q00F_11F = model.vec[251]
    obj$q01F_00F = model.vec[252]
    obj$q01F_10F = model.vec[253]
    obj$q01F_11F = model.vec[254]
    obj$q10F_00F = model.vec[255]
    obj$q10F_01F = model.vec[256]
    obj$q10F_11F = model.vec[257]
    obj$q11F_00F = model.vec[258]
    obj$q11F_01F = model.vec[259]
    obj$q11F_10F = model.vec[260]
    
    obj$q00F_00A = model.vec[261]
    obj$q00F_00B = model.vec[262]
    obj$q00F_00C = model.vec[263]
    obj$q00F_00D = model.vec[264]
    obj$q00F_00E = model.vec[265]
    obj$q00F_00G = model.vec[266]
    obj$q00F_00H = model.vec[267]
    obj$q01F_01A = model.vec[268]
    obj$q01F_01B = model.vec[269]
    obj$q01F_01C = model.vec[270]
    obj$q01F_01D = model.vec[271]
    obj$q01F_01E = model.vec[272]
    obj$q01F_01G = model.vec[273]
    obj$q01F_01H = model.vec[274]
    obj$q10F_10A = model.vec[275]
    obj$q10F_10B = model.vec[276]
    obj$q10F_10C = model.vec[277]
    obj$q10F_10D = model.vec[278]
    obj$q10F_10E = model.vec[279]
    obj$q10F_10G = model.vec[280]
    obj$q10F_10H = model.vec[281]
    obj$q11F_11A = model.vec[282]
    obj$q11F_11B = model.vec[283]
    obj$q11F_11C = model.vec[284]
    obj$q11F_11D = model.vec[285]
    obj$q11F_11E = model.vec[286]
    obj$q11F_11G = model.vec[287]
    obj$q11F_11H = model.vec[288]
    
    ##Hidden State G
    obj$lambda00G = model.vec[289] / (1 + model.vec[293])
    obj$lambda01G = model.vec[290] / (1 + model.vec[294])
    obj$lambda10G = model.vec[291] / (1 + model.vec[295])
    obj$lambda11G = model.vec[292] / (1 + model.vec[296])
    obj$mu00G = (model.vec[293] * model.vec[289]) / (1 + model.vec[293])
    obj$mu01G = (model.vec[294] * model.vec[290]) / (1 + model.vec[294])
    obj$mu10G = (model.vec[295] * model.vec[291]) / (1 + model.vec[295])
    obj$mu11G = (model.vec[296] * model.vec[292]) / (1 + model.vec[296])
    obj$q00G_01G = model.vec[297]
    obj$q00G_10G = model.vec[298]
    obj$q00G_11G = model.vec[299]
    obj$q01G_00G = model.vec[300]
    obj$q01G_10G = model.vec[301]
    obj$q01G_11G = model.vec[302]
    obj$q10G_00G = model.vec[303]
    obj$q10G_01G = model.vec[304]
    obj$q10G_11G = model.vec[305]
    obj$q11G_00G = model.vec[306]
    obj$q11G_01G = model.vec[307]
    obj$q11G_10G = model.vec[308]
    
    obj$q00G_00A = model.vec[309]
    obj$q00G_00B = model.vec[310]
    obj$q00G_00C = model.vec[311]
    obj$q00G_00D = model.vec[312]
    obj$q00G_00E = model.vec[313]
    obj$q00G_00F = model.vec[314]
    obj$q00G_00H = model.vec[315]
    obj$q01G_01A = model.vec[316]
    obj$q01G_01B = model.vec[317]
    obj$q01G_01C = model.vec[318]
    obj$q01G_01D = model.vec[319]
    obj$q01G_01E = model.vec[320]
    obj$q01G_01F = model.vec[321]
    obj$q01G_01H = model.vec[322]
    obj$q10G_10A = model.vec[323]
    obj$q10G_10B = model.vec[324]
    obj$q10G_10C = model.vec[325]
    obj$q10G_10D = model.vec[326]
    obj$q10G_10E = model.vec[327]
    obj$q10G_10F = model.vec[328]
    obj$q10G_10H = model.vec[329]
    obj$q11G_11A = model.vec[330]
    obj$q11G_11B = model.vec[331]
    obj$q11G_11C = model.vec[332]
    obj$q11G_11D = model.vec[333]
    obj$q11G_11E = model.vec[334]
    obj$q11G_11F = model.vec[335]
    obj$q11G_11H = model.vec[336]
    
    ##Hidden State H
    obj$lambda00H = model.vec[337] / (1 + model.vec[341])
    obj$lambda01H = model.vec[338] / (1 + model.vec[342])
    obj$lambda10H = model.vec[339] / (1 + model.vec[343])
    obj$lambda11H = model.vec[340] / (1 + model.vec[344])
    obj$mu00H = (model.vec[341] * model.vec[337]) / (1 + model.vec[341])
    obj$mu01H = (model.vec[342] * model.vec[338]) / (1 + model.vec[342])
    obj$mu10H = (model.vec[343] * model.vec[339]) / (1 + model.vec[343])
    obj$mu11H = (model.vec[344] * model.vec[340]) / (1 + model.vec[344])
    obj$q00H_01H = model.vec[345]
    obj$q00H_10H = model.vec[346]
    obj$q00H_11H = model.vec[347]
    obj$q01H_00H = model.vec[348]
    obj$q01H_10H = model.vec[349]
    obj$q01H_11H = model.vec[350]
    obj$q10H_00H = model.vec[351]
    obj$q10H_01H = model.vec[352]
    obj$q10H_11H = model.vec[353]
    obj$q11H_00H = model.vec[354]
    obj$q11H_01H = model.vec[355]
    obj$q11H_10H = model.vec[356]
    
    obj$q00H_00A = model.vec[357]
    obj$q00H_00B = model.vec[358]
    obj$q00H_00C = model.vec[359]
    obj$q00H_00D = model.vec[360]
    obj$q00H_00E = model.vec[361]
    obj$q00H_00F = model.vec[362]
    obj$q00H_00G = model.vec[363]
    obj$q01H_01A = model.vec[364]
    obj$q01H_01B = model.vec[365]
    obj$q01H_01C = model.vec[366]
    obj$q01H_01D = model.vec[367]
    obj$q01H_01E = model.vec[368]
    obj$q01H_01F = model.vec[369]
    obj$q01H_01G = model.vec[370]
    obj$q10H_10A = model.vec[371]
    obj$q10H_10B = model.vec[372]
    obj$q10H_10C = model.vec[373]
    obj$q10H_10D = model.vec[374]
    obj$q10H_10E = model.vec[375]
    obj$q10H_10F = model.vec[376]
    obj$q10H_10G = model.vec[377]
    obj$q11H_11A = model.vec[378]
    obj$q11H_11B = model.vec[379]
    obj$q11H_11C = model.vec[380]
    obj$q11H_11D = model.vec[381]
    obj$q11H_11E = model.vec[382]
    obj$q11H_11F = model.vec[383]
    obj$q11H_11G = model.vec[384]
    
    obj$psi <- model.vec[385]
    
    return(obj)
}


print.muhisse.fit <- function(x,...){
    ## Function to print a "muhisse.fit" object.
    set.zero <- max( x$index.par )
    ## Keep only the parameters estimated:
    par.list <- x$solution[ !x$index.par == set.zero ]
    ntips <- Ntip( x$phy )
    nstates <- ncol( x$trans.matrix )/4
    output <- c(x$loglik, x$AIC, x$AICc, ntips, nstates)
    names(output) <- c("lnL", "AIC", "AICc", "n.taxa", "n.hidden.states")
    cat("\n")
    cat("Fit \n")
    print(output)
    cat("\n")
    cat("Model parameters: \n")
    cat("\n")
    print(par.list)
    cat("\n")
    if(!is.null(x$psi.type)){
        if(x$psi.type == "m_only"){
            cat("psi estimate reflects fossil tip sampling only. \n")
        }else{
            cat("psi estimate reflects both fossil edge and tip sampling. \n")
        }
    }
}

Try the hisse package in your browser

Any scripts or data that you put into this service are public.

hisse documentation built on Feb. 16, 2023, 10:26 p.m.