gn3sim<-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)
    height <- c(theta[17:length(theta)], 1)
    merge <- merge2
    k <- dim(merge)[1]
    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 <- merge[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 <- genseq2(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])
    }
}

genseq2 <- function(t1, t2, seq, Sx2, Sy2, Pix, Piy) {
    sequ <- matrix(0, nc = 2, nr = length(seq))
    p1 <- Pt(Sx2, Pix, t1)
    p2 <- Pt(Sy2, Piy, t2)
    for(i in 1:length(seq)) {
        sequ[i, 1] <- sample(1:4, 1, prob = p1[seq[i],  ])
        sequ[i, 2] <- sample(1:4, 1, prob = p2[seq[i],  ])
    }
    return(sequ)
}