simemb<-function(theta, seqLength, merge2)
{
    pix1 <- theta[1:3]
    piy1 <- theta[4:6]
    f00 <- theta[7:9]
    ss <- theta[10:15]
    rho <- theta[16]
    pix2 <- c(pix1, 1 - sum(pix1))
    Pix2 <- diag(pix2)
    piy2 <- c(piy1, 1 - sum(piy1))
    Piy2 <- diag(piy2)
    f0 <- c(f00, 1 - sum(f00))
    Sx1 <- Smatrix(ss, pix2)
    Sy1 <- rho * Smatrix(ss, piy2)
    merge <- merge2
    k <- dim(merge)[1]
    height2 <- c(theta[17:length(theta)], 1)
    height <- height2
    h <- c(0, height)
    me <- merge[k,  ]
    m <- me
    t0 <- h[k + 1]
    cc <- 0
    if(merge[, 1][k] > 0) {
        cc <- merge[, 1][k]
        v <- matrix(0, k, 1)
        v[1] <- merge[k, 1]
        cur <- 0
        max <- 1
        while(cur < max) {
            cur <- cur + 1
            for(i in c(1, 2)) {
                if(merge[(v[cur]), i] > 0) {
                    max <- max + 1
                    v[max] <- merge[(v[cur]),
                        i]
                }
            }
            cc <- matrix(0, max, 1)
            for(i in 1:max) {
                cc[i] <- v[i]
            }
        }
    }
    seqr <- matrix(0, seqLength, k - 1)
    fs <- matrix(0, seqLength, k + 1)
    seq1 <- cbind(sample(1:4, seqLength, prob = f0, replace = T)
        )
    for(i in 1:k) {
        iii <- k - i + 1
        if(iii == k) {
            xx <- 0
            if(any(merge[k, 1] == cc)) {
                pix <- pix2
                Pix <- Pix2
                piy <- piy2
                Piy <- Piy2
                Sx2 <- Sx1
                Sy2 <- Sy1
            }
            else {
                piy <- pix2
                Piy <- Pix2
                pix <- piy2
                Pix <- Piy2
                Sy2 <- Sx1
                Sx2 <- Sy1
            }
        }
        else {
            if(any(iii == cc)) {
                pix <- pix2
                Pix <- Pix2
                piy <- pix2
                Piy <- Pix2
                Sx2 <- Sx1
                Sy2 <- Sx1
            }
            else {
                piy <- piy2
                Piy <- Piy2
                pix <- piy2
                Pix <- Piy2
                Sx2 <- Sy1
                Sy2 <- Sy1
            }
        }
        m <- rev(sort(me))
        mm <- merge2[m[1],  ]
        if(me[1] > 0) {
            t1 <- t0 - h[me[1] + 1]
        }
        else {
            t1 <- t0
        }
        if(me[2] > 0) {
            t2 <- t0 - h[me[2] + 1]
        }
        else {
            t2 <- t0
        }
        if(t1 < 0 || t1 > 1) {
            return(0)
        }
        if(t2 < 0 || t2 > 1) {
            return(0)
        }
        hp <- rev(order(me))
        if(i != 1) {
            seq1 <- seqr[, (k - i + 1)]
        }
        seq2 <- genseq4(t1, t2, seq1, Sx2, Sy2, Pix, Piy)
        if(me[1] < 0) {
            fs[, (abs(me[1]))] <- seq2[, 1]
        }
        else {
            seqr[, (abs(me[1]))] <- seq2[, 1]
        }
        if(me[2] < 0) {
            fs[, (abs(me[2]))] <- seq2[, 2]
        }
        else {
            seqr[, (abs(me[2]))] <- seq2[, 2]
        }
        if(m[1] > 0)
            t0 <- h[m[1] + 1]
        else {
            return(fs)
        }
        me <- c(mm, m[-1])
    }
}

genseq4 <- function(t1, t2, seq, Sx2, Sy2, Pix, Piy) {
    seqx <- matrix(0, nc = 1, nr = length(seq))
    seqy <- matrix(0, nc = 1, nr = length(seq))
    Rx <- Sx2 %*% Pix
    Ry <- Sy2 %*% Piy
    for(j in 1:(length(seq))) {
        Tx <- matrix(0, 1, 10)
        WTx <- matrix(0, 1, 10)
        m <- matrix(0, 1, 10)
        m[1] <- seq[j]
        WTx[1] <- rexp(1, ( - Rx[m[1], m[1]]))
        xx <- m[1]
        for(i in 2:10) {
            Tx[i] <- Tx[i - 1] + WTx[i - 1]
            if(Tx[i] >= t1)
                break
            rrx <- Rx[m[i - 1],  ]/( - Rx[m[i - 1], m[i - 1]])
            rrx[m[i - 1]] <- 0
            m[i] <- sample(1:4, 1, rrx, replace = T)
            WTx[i] <- rexp(1, ( - Rx[m[i], m[i]]))
            xx <- m[i]
        }
        seqx[j] <- xx
    }
    for(k in 1:(length(seq))) {
        Ty <- matrix(0, 1, 10)
        WTy <- matrix(0, 1, 10)
        n <- matrix(0, 1, 10)
        n[1] <- seq[k]
        WTy[1] <- rexp(1, ( - Rx[n[1], n[1]]))
        yy <- n[1]
        for(i in 2:10) {
            Ty[i] <- Ty[i - 1] + WTy[i - 1]
            if(Ty[i] >= t2)
                break
            rry <- Ry[n[i - 1],  ]/( - Ry[n[i - 1], n[i - 1]])
            rry[n[i - 1]] <- 0
            n[i] <- sample(1:4, 1, rry, replace = T)
            WTy[i] <- rexp(1, ( - Ry[n[i], n[i]]))
            yy <- n[i]
        }
        seqy[k] <- yy
    }
    sequ <- cbind(seqx, seqy)
    return(sequ)
}