gn<-function(theta, 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)) height <- theta[17:length(theta)] if(sum(abs(pix2)) > 1) { return(0) } if(sum(abs(piy2)) > 1) { return(0) } if(sum(abs(f0)) > 1) { return(0) } if(rho <= 0) { return(0) } if(any(theta <= 0)) { return(0) } Sx1 <- Smatrix(ss, pix2) Sy1 <- rho * Smatrix(ss, piy2) k <- dim(merge2)[1] if(k == 1) { F1 <- Fmatrix(theta[17], theta[17], f0, Sx1, Sy1, Pix2, Piy2) return(F1) } else { h <- c(0, height, 1) me <- merge2[k, ] f1 <- f0 m <- me t0 <- h[k + 1] cc <- 0 if(merge2[, 1][k] > 0) { cc <- merge2[, 1][k] v <- matrix(0, k, 1) v[1] <- merge2[k, 1] cur <- 0 max <- 1 while(cur < max) { cur <- cur + 1 for(i in c(1, 2)) { if(merge2[(v[cur]), i] > 0) { max <- max + 1 v[max] <- merge2[ (v[cur]), i] } } cc <- matrix(0, max, 1) for(i in 1:max) { cc[i] <- v[i] } } } for(i in 1:k) { zzz <- 0 iii <- k - i + 1 if(iii == k) { xx <- 0 if(any(merge2[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)) for(j in 1:(4^(i - 1))) { if(j == 1) { F1 <- Fmatrix(t1, t2, f1[ 1:4], Sx2, Sy2, Pix, Piy) } else { F1 <- Fmatrix(t1, t2, f1[ (4 * (j - 1) + 1): (4 * (j - 1) + 4)], Sx2, Sy2, Pix, Piy) } if(j == 1) { f <- c(as.vector(F1), f1[ -1:-4]) } else { f <- c(f[1:(4^2 * (j - 1))], as.vector(F1), f1[ -1: - (4 * j)]) } } f12 <- array(f, c(rep(4, i + 1))) f112 <- aperm(f12, c(hp)) f1 <- as.vector(f112) if(m[1] > 0) { t0 <- h[m[1] + 1] } else { return(f112) } me <- c(mm, m[-1]) } } }