#   Two-parameter likelihood plots

loglik.plot <- function(panel) {
   panel$par1.low  <- as.numeric(panel$par1.low)
   panel$par1.high <- as.numeric(panel$par1.high)
   panel$par2.low  <- as.numeric(panel$par2.low)
   panel$par2.high <- as.numeric(panel$par2.high)
      par1       <- seq(panel$par1.low, panel$par1.high, length = panel$ngrid)
      par2       <- seq(panel$par2.low, panel$par2.high, length = panel$ngrid)
      parmat     <- as.matrix(expand.grid(par1, par2))
      loglik.mat <- apply(parmat, 1, panel$loglik.fn, data = panel$data)
      loglik.mat <- matrix(loglik.mat, ncol = panel$ngrid)
      if (panel$transparent) alpha <- rep(0.5, panel$ngrid^2)
      else                   alpha <- rep(  1, panel$ngrid^2)
      if (panel$lose.lower.region) {
         ind <- (loglik.mat <= max(loglik.mat) - 10)
         loglik.mat[ind] <- max(loglik.mat) - 10
         alpha[c(ind)] <- 0
         }
      ints <- seq(min(loglik.mat), max(loglik.mat), length = 51)
      ind  <- findInterval(c(loglik.mat), ints)
      clr  <- terrain.colors(50)[ind]
      scaling <- rp.plot3d.new(range(par1), range(loglik.mat), range(par2), 
                     xlab = "par 1", ylab = "log-likelihood", zlab = "par 2",
                     type = "n", new.window = FALSE)
      ax <- scaling(par1,  par1,  par1)$x
      ay <- scaling(loglik.mat, loglik.mat, loglik.mat)$y
      az <- scaling(par2, par2, par2)$z
      rgl.surface(ax, az, ay, alpha = alpha, col = clr)
      material3d(alpha = 1)
      panel$loglik.mat <- loglik.mat
      panel$scaling <- scaling
      panel$n.add <- 0
   panel <- loglik.plot.add(panel)
   panel
   }

loglik.plot.add <- function(panel) {
   panel$par1.low  <- as.numeric(panel$par1.low)
   panel$par1.high <- as.numeric(panel$par1.high)
   panel$par2.low  <- as.numeric(panel$par2.low)
   panel$par2.high <- as.numeric(panel$par2.high)
   if (panel$n.add > 0) {
      for (i in 1:panel$n.add) rgl.pop()
      panel$n.add <- 0
      }
   with(panel, {
      if (mle.showing) {
         if ((mloglik$par[1] >= par1.low) & (mloglik$par[1] <= par1.high) & 
             (mloglik$par[2] >= par2.low) & (mloglik$par[2] <= par2.high)) {
            a <- scaling(mloglik$par[1], mloglik$value, mloglik$par[2])
            points3d(a$x, a$y, a$z, col = "red", size = 5)
            }
         }
      if (ci.showing) {
         a <- scaling(mloglik$par[1], mloglik$value, mloglik$par[2])
         if ((mloglik$value - 3 > min(loglik.mat)) & 
             (mloglik$value - 3 < max(loglik.mat))) {
            a <- scaling(c(par1.low, par1.low, par1.high, par1.high),
                   rep(mloglik$value - 3, 4), 
                   c(par2.high, par2.low, par2.low, par2.high))
            quads3d(a$x, a$y, a$z, col = "red", size = 5)
            }
         }
      if (quad.showing & ("hessian" %in% names(mloglik))) {
         par1    <- seq(par1.low, par1.high, length = ngrid)
         par2    <- seq(par2.low, par2.high, length = ngrid)
         parmat  <- as.matrix(expand.grid(par1, par2))
         quad.fn <- function(par, mloglik) mloglik$value + 
              0.5 * as.vector((par - mloglik$par) %*% mloglik$hessian %*% (par - mloglik$par))
         loglik.mat    <- apply(parmat, 1, quad.fn, mloglik = mloglik)
         loglik.mat    <- matrix(loglik.mat, ncol = ngrid)
         if (transparent) alpha <- rep(0.5, ngrid^2)
            else          alpha <- rep(  1, ngrid^2)
         ind <- (loglik.mat <= max(loglik.mat) - 10)
         loglik.mat[ind] <- max(loglik.mat) - 10
         alpha[c(ind)] <- 0
         ints <- seq(min(loglik.mat), max(loglik.mat), length = 51)
         ind  <- findInterval(c(loglik.mat), ints)
         clr  <- terrain.colors(50)[ind]
         ax <- scaling(par1,  par1,  par1)$x
         ay <- scaling(loglik.mat, loglik.mat, loglik.mat)$y
         az <- scaling(par2, par2, par2)$z
         rgl.surface(ax, az, ay, alpha = alpha, col = clr, front = "line", back = "line")
         material3d(alpha = 1)
         }
      })
   panel$n.add <- 0
   if (panel$mle.showing & (panel$mloglik$par[1] >= panel$par1.low)  & 
                           (panel$mloglik$par[1] <= panel$par1.high) & 
                           (panel$mloglik$par[2] >= panel$par2.low)  & 
                           (panel$mloglik$par[2] <= panel$par2.high))
      panel$n.add <- panel$n.add + 1
   if (panel$ci.showing & (panel$mloglik$value - 3 > min(panel$loglik.mat)) & 
                          (panel$mloglik$value - 3 < max(panel$loglik.mat)))   
      panel$n.add <- panel$n.add + 1
   if (panel$quad.showing & ("hessian" %in% names(panel$mloglik))) 
      panel$n.add <- panel$n.add + 1
   panel
   }

help.fn <- function(panel) {
   system("acroread temp.pdf &")
   panel
   }

rp.loglik2 <- function(par1.range, par2.range, loglik.fn, data) {
   rgl.open()
   rgl.bg(col = c("white", "black"))
   mloglik <- optim(c(mean(par1.range), mean(par2.range)), loglik.fn, data = data,
                       control = list(fnscale = -1), hessian = TRUE)
   # print(mloglik$par)
   # print(mloglik$value)
   if (mloglik$convergence > 0)
      print(mloglik)
   panel <- rp.control(data = data, ngrid = 50, n.add = 0, mloglik = mloglik,
                       loglik.fn = loglik.fn,
                       par1.low  = par1.range[1], 
                       par1.high = par1.range[2],
                       par2.low  = par2.range[1], 
                       par2.high = par2.range[2],
                       lose.lower.region = TRUE)
   rp.textentry(panel, par1.low,          loglik.plot)
   rp.textentry(panel, par1.high,         loglik.plot)
   rp.textentry(panel, par2.low,          loglik.plot)
   rp.textentry(panel, par2.high,         loglik.plot)
   rp.checkbox(panel,  mle.showing,       loglik.plot.add)
   rp.checkbox(panel,  ci.showing,        loglik.plot.add)
   rp.checkbox(panel,  quad.showing,      loglik.plot.add)
   rp.checkbox(panel,  lose.lower.region, loglik.plot)
   rp.checkbox(panel,  transparent,       loglik.plot)
   # rp.button(panel, help.fn, title = "Help")
   rp.do(panel, loglik.plot)
   }

rp.plot3d.new <- function (x, y, z, xlab = NA, ylab = NA, zlab =
NA, axes = TRUE, new.window = TRUE,
    type = "p", size = 3, col = "red", xlim = NA, ylim = NA,
    zlim = NA, ...)
{
    if (require(rgl)) {
        if (is.na(xlab))
            xlab <- deparse(substitute(x))
        if (is.na(ylab))
            ylab <- deparse(substitute(y))
        if (is.na(zlab))
            zlab <- deparse(substitute(z))
        xrange <- xlim
        yrange <- ylim
        zrange <- zlim
        ind <- !is.na(x + y + z)
        if (length(col) == length(x))
            ind <- (ind & (!is.na(col)))
        if (!all(ind))
            cat("Warning: missing data removed. \n")
        if (any(is.na(xlim))) {
            xrange[1] <- min(x[ind]) - 0.05 * diff(range(x[ind]))
            xrange[2] <- max(x[ind]) + 0.05 * diff(range(x[ind]))
        }
        if (any(is.na(ylim))) {
            yrange[1] <- min(y[ind]) - 0.05 * diff(range(y[ind]))
            yrange[2] <- max(y[ind]) + 0.05 * diff(range(y[ind]))
        }
        if (any(is.na(zlim))) {
            zrange[1] <- min(z[ind]) - 0.05 * diff(range(z[ind]))
            zrange[2] <- max(z[ind]) + 0.05 * diff(range(z[ind]))
        }
        xscale <- pretty(xrange)
        yscale <- pretty(yrange)
        zscale <- pretty(zrange)
        xscale <- xscale[xscale >= xrange[1] & xscale <= xrange[2]]
        yscale <- yscale[yscale >= yrange[1] & yscale <= yrange[2]]
        zscale <- zscale[zscale >= zrange[1] & zscale <= zrange[2]]
        xadj1 <- mean(xrange)
        yadj1 <- mean(yrange)
        zadj1 <- mean(zrange)
        xadj2 <- diff(xrange)/2
        yadj2 <- diff(yrange)/2
        zadj2 <- diff(zrange)/2
        x.orig <- x
        y.orig <- y
        z.orig <- z
        x <- (x - xadj1)/xadj2
        y <- (y - yadj1)/yadj2
        z <- (z - zadj1)/zadj2
        xscale.adj <- (xscale - xadj1)/xadj2
        yscale.adj <- (yscale - yadj1)/yadj2
        zscale.adj <- (zscale - zadj1)/zadj2
        rx <- c(-1, 1)
        ry <- c(-1, 1)
        rz <- c(-1, 1)
        if (new.window) {
           rgl.open()
           rgl.bg(col = c("white", "black"))
           }
        else
           rgl.clear()
        rgl.viewpoint(-30, 30, fov = 1)
        if (axes) {
            rgl.lines(rx[c(1, 2, 2, 2, 2, 1, 1, 1)], ry[rep(1,
                8)], rz[c(1, 1, 1, 2, 2, 2, 2, 1)], col = "black")
            rgl.lines(rx[c(1, 2, 2, 2, 2, 1, 1, 1)], ry[rep(2,
                8)], rz[c(1, 1, 1, 2, 2, 2, 2, 1)], col = "black")
            for (i in 1:2) for (j in 1:2) rgl.lines(rx[c(i, i)],
                ry[c(1, 2)], rz[c(j, j)], col = "black")
            rgl.texts(mean(rx), min(rx), min(rx), "")
            delta <- 0.1
            nyticks <- length(yscale)
            if (nyticks/2 - floor(nyticks/2) > 0)
                ypos <- 1/(nyticks - 1)
            else ypos <- 0
            rgl.texts(c(0, -1 - 2 * delta, -1 - 2 * delta), c(-1 -
                2 * delta, ypos, -1 - 2 * delta), c(1 + 2 * delta,
                -1 - 2 * delta, 0), c(xlab, ylab, zlab), adj = c(0.5,
                0.5, 1), col = "blue")
            rgl.texts((xscale - xadj1)/xadj2, -1 - delta, 1 +
                delta, as.character(xscale), col = "black")
            rgl.texts(-1 - delta, (yscale - yadj1)/yadj2, -1 -
                delta, as.character(yscale), col = "black")
            rgl.texts(-1 - delta, -1 - delta, (zscale - zadj1)/zadj2,
                as.character(zscale), col = "black")
            scaling <- function(x, y, z) {
                list(x = x, y = y, z = z)
            }
            rgl.segments(xscale.adj, -1, 1, xscale.adj, -1 -
                delta/4, 1 + delta/4, scaling = scaling, col = "black")
            rgl.segments(-1, yscale.adj, -1, -1 - delta/4, yscale.adj,
                -1 - delta/4, scaling = scaling, col = "black")
            rgl.segments(-1, -1, zscale.adj, -1 - delta/4, -1 -
                delta/4, zscale.adj, scaling = scaling, col = "black")
            }
        if (!(type == "n")) {
           ind1 <- ((x.orig >= xrange[1]) & (x.orig <= xrange[2]) &
               (y.orig >= yrange[1]) & (y.orig <= yrange[2]) &
               (z.orig >= zrange[1]) & (z.orig <= zrange[2]))
           ind <- (ind1 & ind)
           if (length(col) == length(x.orig))
               clr <- col[ind]
           else clr <- col
           rgl.points(x[ind], y[ind], z[ind], size = size, col = clr)
           }
        scaling <- function(x, y, z) {
            xx <- (x - xadj1)/xadj2
            yy <- (y - yadj1)/yadj2
            zz <- (z - zadj1)/zadj2
            list(x = xx, y = yy, z = zz)
        }
        invisible(scaling)
    }
    else {
        warning("Package RGL is not installed.")
    }
}

rgl.segments <- function(x0, y0, z0, x1, y1, z1, scaling, ...) {
         a <- scaling(c(rbind(x0, x1)), c(rbind(y0, y1)), c(rbind(z0, z1)))
         rgl.lines(a$x, a$y, a$z, ...)
         }