"ellipse" <-
  function (x, ...) 
  UseMethod("ellipse")
"ellipse.default" <-
  function (x, scale = c(1, 1), centre = c(0, 0), level = 0.95, 
            t = sqrt(qchisq(level, 2)), which = c(1, 2), npoints = 100) 
{
  names <- c("x", "y")
  if (is.matrix(x)) {
    xind <- which[1]
    yind <- which[2]
    scale <- sqrt(c(x[xind, xind], x[yind, yind]))
    r <- x[xind, yind]/scale[1]/scale[2]
    if (!is.null(dimnames(x)[[1]])) 
      names <- dimnames(x)[[1]][c(xind, yind)]
  }
  else r <- x
  d <- acos(r)
  a <- seq(0, 2 * pi, len = npoints)
  matrix(c(t * scale[1] * cos(a + d/2) + centre[1], t * scale[2] * 
           cos(a - d/2) + centre[2]), npoints, 2, dimnames = list(NULL, 
                                                    names))
}
"ellipse.glm" <-
  function (x, which = c(1, 2), level = 0.95, t, ...) 
{
  s <- summary(x)
  if (missing(t)) {
    t <- switch(as.character(x$family["family"]), binomial = sqrt(qchisq(level, 2)), poisson = sqrt(qchisq(level, 2)), sqrt(2 * qf(level, 2, s$df[2])))
  }
  ellipse.default(s$dispersion * s$cov.unscaled[which, which], 
                  centre = x$coefficients[which], t = t)
}
"ellipse.lm" <-
  function (x, which = c(1, 2), level = 0.95, t = sqrt(2 * qf(level, 
                                                2, x$df.residual)), ...) 
{
  s <- summary(x)
  ellipse.default(s$sigma^2 * s$cov.unscaled[which, which], 
                  centre = x$coefficients[which], t = t)
}
"ellipse.nls" <-
  function (x, which = c(1, 2), level = 0.95, t = sqrt(2 * qf(level, 
                                                2, s$df[2])), ...) 
{
  s <- summary(x)
  ellipse.default(s$sigma^2 * s$cov.unscaled[which, which], 
                  centre = x$m$getPars()[which], t = t)
}
"ellipse.profile" <-
  function (x, which = c(1, 2), level = 0.95, t = sqrt(qchisq(level, 
                                                2)), npoints = 100) 
{
  aa <- x[[which[1]]][[2]][, which[1]]
  ar <- x[[which[1]]][[2]][, which[2]]
  ra <- x[[which[2]]][[2]][, which[1]]
  rr <- x[[which[2]]][[2]][, which[2]]
  atau <- x[[which[1]]][[1]]
  rtau <- x[[which[2]]][[1]]
  arange <- range(c(aa, ra))
  rrange <- range(c(ar, rr))
  atau <- atau/t
  rtau <- rtau/t
  getad <- function(tau1, tau2) {
    if (abs(tau1) > 1) 
      tau1 <- tau1/abs(tau1)
    if (abs(tau2) > 1) 
      tau2 <- tau2/abs(tau2)
    acos1 <- acos(tau1)
    acos2 <- acos(tau2)
    d <- abs(acos1 - acos2)
    a <- (acos1 + acos2)/2
    if (acos1 < acos2) 
      a <- -a
    c(a, d)
  }
  myapprox <- function(x, y, where) {
    good <- is.finite(x) & is.finite(y)
    x <- x[good]
    y <- y[good]
    if (length(x) > 1) {
      result <- approx(x[good], y[good], where)$y
      bad <- is.na(result)
      if (any(bad)) {
        for (i in 1:length(result)) {
          if (bad[i]) {
            if (where[i] > x[length(x)]) {
              x1 <- x[length(x) - 1]
              y1 <- y[length(x) - 1]
              x2 <- x[length(x)]
              y2 <- y[length(x)]
            }
            else if (where[i] < x[1]) {
              x1 <- x[1]
              y1 <- y[1]
              x2 <- x[2]
              y2 <- y[2]
            }
            else stop("Unexpected NA")
            result[i] <- y1 + (where[i] - x1)/(x2 - x1) * 
              (y2 - y1)
          }
        }
      }
    }
    else result <- rep(y, length(where))
    result
  }
  ad <- matrix(NA, nrow = 5, ncol = 2)
  ad[1, ] <- getad(1, myapprox(rr, rtau, myapprox(aa, ar, myapprox(atau,aa, 1))))
  ad[2, ] <- getad(myapprox(aa, atau, myapprox(rr, ra, myapprox(rtau, rr, 1))), 1)
  ad[3, ] <- getad(-1, myapprox(rr, rtau, myapprox(aa, ar, myapprox(atau, aa, -1))))
  ad[4, ] <- getad(myapprox(aa, atau, myapprox(rr, ra, myapprox(rtau, rr, -1))), -1)
  i <- order(ad[1:4, 1])
  ad[1:4, ] <- ad[i, ]
  ad[5, 1] <- ad[1, 1] + 2 * pi
  ad[5, 2] <- ad[1, 2]
  ad <- ad[!duplicated(ad[, 1]), ]
  adv <- spline(ad,n = npoints,method= "periodic")
  avals <- adv$x
  dvals <- adv$y
  matrix(c(myapprox(atau, aa, cos(avals + dvals/2)), myapprox(rtau,  rr, cos(avals - dvals/2))), length(avals), 2, dimnames = list(NULL, names(x[which])))
}
"ellipse.profile.glm" <-
  function (x, level = 0.95, t, ...) 
{
  if (missing(t)) {
    t <- switch(attr(x, "original.fit")$family$family, binomial = sqrt(qchisq(level,2)), poisson = sqrt(qchisq(level, 2)), sqrt(2 * qf(level,2, attr(attr(x,"original.fit"), "summary")$df[2])))
  }
  ellipse.profile(x, level = level, t = t, ...)
}

"ellipse.profile.nls" <-
  function (x, level = 0.95, t = sqrt(2 * qf(level, 2, attr(x, 
                               "summary")$df[2])), ...) 
{
  ellipse.profile(x, level = level, t = t, ...)
}
"pairs.profile.nls" <-
  function (x, labels = c(names(x), "Profile tau"), panel = lines, 
            invert = TRUE, plot.tau = T, plot.trace = T, plot.sketch = T, 
            plot.ellipse = F, level = 0.95, ...) 
{
  doaxis <- function(which, dolabel = T) axis(which,labels=dolabel) # outer = T, line = -0.5, labels = dolabel)
  setup <- function(x, y, ...) plot(range(x[!is.na(x)]), 
                                       range(y[!is.na(y)]), type = "n", axes = F, ...)
  if (is.character(panel)) 
    panel <- get(panel, mode = "function")
  n <- length(x)
  if (plot.tau) 
    n <- n + 1
  oldpar <- par("oma", "mar", "cex", "tck", "mgp", "mex", 
                "mfrow")
  oldcex <- par("cex")
  CEX <- oldcex * max(7.7/(2 * n + 3), 0.6)
  par(mfrow = c(n, n), mgp = c(2, 0.8, 0), oma = rep(3, 4), 
      mar = rep(0.5, 4), tck = -0.03/n)
  on.exit({
    par(oldpar)
  })
  par(cex = CEX)
  if (length(labels) < n) 
    labels <- paste(deparse(substitute(x)), "[,", 1:n, "]", 
                    sep = "")
  if (par("pty") == "s") {
    dif <- diff(par("fin"))/2
    if (dif > 0) 
      par(omi = c(dif * n, 0, dif * n, 0) + par("omi"))
    else par(omi = c(0, (-dif) * n, 0, (-dif) * n) + par("omi"))
  }
  alltau <- unlist(lapply(x, function(x) x[[1]]), use.names = F)
  order <- if (invert) 
    1:n
  else n:1
  for (i in order) {
    for (j in 1:n) {
      if (i<=length(x))
          icomp <- x[[i]]
      if (j<=length(x))
          jcomp <- x[[j]]
      xx1 <- NA
      xx2 <- NA
      yy1 <- NA
      yy2 <- NA
      if (i <= length(x)) {
        yy1 <- icomp[[2]][, i]
        if (j <= length(x)) {
          xx1 <- icomp[[2]][, j]
          xx2 <- jcomp[[2]][, j]
          yy2 <- jcomp[[2]][, i]
        }
        else {
          xx1 <- icomp[[1]]
        }
      }
      else {
        yy1 <- jcomp[[1]]
        if (j <= length(x)) {
          xx1 <- jcomp[[2]][, j]
        }
      }
      xx <- c(xx1, NA, xx2)
      yy <- c(yy1, NA, yy2)
      if (i <= length(x)) {
        if (j <= length(x)) 
          setup(xx, yy, ...)
        else setup(alltau, yy, ...)
      }
      else {
        if (j <= length(x)) 
          setup(xx, alltau, ...)
        else setup(alltau, alltau)
      }
      box()
      if (i == 1) 
        doaxis(3, j%%2 == 0)
      if (i == n) 
        doaxis(1, j%%2 == 1)
      if (j == 1) 
        doaxis(2, i%%2 == 0)
      if (j == n) 
        doaxis(4, i%%2 == 1)
      if (i != j) {
        if ((i <= length(x)) && (j <= length(x))) {
          if (plot.trace) 
            panel(xx, yy, ...)
          if (plot.sketch) 
            for (l in level) panel(ellipse(x, which = c(j, 
                                                i), level = l), ...)
          if (plot.ellipse && !is.null(fit <- attr(x, 
                                                   "original.fit"))) 
            for (l in level) panel(ellipse(fit, which = c(j, 
                                                  i), level = l), ...)
        }
        else if (plot.tau) 
          panel(xx, yy, ...)
      }
      else {
        par(usr = c(0, 1, 0, 1))
        text(0.5, 0.5, labels[i], cex = 1.5 * CEX)
      }
    }
  }
  invisible()
}
"plotcorr" <-
  function (corr, outline = T, dev = T, col = 4, paropts = NULL, 
            numbers = F, ...) 
{
  if (deparse(substitute(dev)) == "postscript") {
    cat("feature not still implemented in ellipse for R, sorry\n")
#    postscript(preamble = ps.preamble.ellipse, font = ps.fonts.ellipse, 
#               ...)
#    assign("ellipse.fontnum", length(ps.fonts.ellipse), where = 0)
#    cat("postscript device started with ellipse fonts\n")
  }
  else if (is.function(dev)) 
    dev(...)
  par(pty = "s", mar = c(5, 0, 4, 0) + 0.1)
#  par(pty = "s", mar = c(5, 4, 4, 2) + 0.1)

  if (!is.null(paropts)) 
    par(paropts)
  if (is.null(corr)) 
    return(invisible())
  if ((!is.matrix(corr)) || (round(min(corr), 6) < -1) || (round(max(corr), 
                                                                 6) > 1)) 
    stop("Need a correlation matrix")
  rowdim <- dim(corr)[1]
  coldim <- dim(corr)[2]
  maxdim <- max(rowdim, coldim)
  rowlabs <- dimnames(corr)[[1]]
  collabs <- dimnames(corr)[[2]]
  if (is.null(rowlabs)) 
    rowlabs <- 1:rowdim
  if (is.null(collabs)) 
    collabs <- 1:coldim
#1st change
#cxy <- par("cxy")
#cxy[1]<-yinch(par("csi"))
#cxy[2]<-yinch(par("csi"))
#cxy<-c(0.57*xinch(par("csi")), yinch(par("csi")))
  cxy<-par("cin")/par("pin")
  xlabwidth <- (max(nchar(rowlabs)) * cxy[1])/cxy[2]
  ylabwidth <- (max(nchar(collabs)) * cxy[1])/cxy[2]
#  plot(c(-xlabwidth, maxdim + 0.5), c(0.5, maxdim + 1 + ylabwidth), 
#       type = "n", bty = "n", axes = F, xlab = "", ylab = "",asp=1)
  plot(c(-xlabwidth, maxdim + 0.5), c(0.5, maxdim + 1 + ylabwidth), 
       type = "n", bty = "n", axes = F, xlab = "", ylab = "",asp=1)
                                        #2nd change
                                        #par(cex=0.57*par("cex"))
                                        #par(cex = par("cex")/(par("cxy")[2]))
  text(rep(0, rowdim), rowdim:1, labels = rowlabs, adj = 1)
  text(1:coldim, rep(rowdim + 1, coldim), labels = collabs, 
       srt = 90, adj = 0)
  cols <- rep(1:coldim, rep(rowdim, coldim))
  rows <- rep(1:rowdim, coldim)
  if (!numbers) {
#    if (names(dev.cur()) == "postscript") {
#      text(cols, rowdim + 1 - rows, ellipse.chars[round((as.vector(corr) + 
#                                                         1)/2 * (length(ellipse.chars) - 1) + 1)], 
#           col = col)
#      if (outline) 
#        text(cols, rowdim + 1 - rows, ellipse.chars[round((as.vector(corr) + 
#                                                           1)/2 * (length(ellipse.chars) - 1) + 1)])
#    }
#    else {
      mat <- diag(c(1, 1))
      for (i in 1:length(cols)) {
        mat[1, 2] <- as.vector(corr)[i]
        mat[2, 1] <- mat[1, 2]
        ell <- ellipse(mat, t = 0.43)
        ell[, 1] <- ell[, 1] + cols[i]
        ell[, 2] <- ell[, 2] + rowdim + 1 - rows[i]
        polygon(ell, col = col)
        if (outline) 
          lines(ell)
      }
#    }
  }
  else text(cols + 0.3, rowdim + 1 - rows, round(10 * as.vector(corr), 
                                                 0), adj = 1, cex = 0.75 * par("cex"))
  invisible()
}







"profile.glm" <-
  function (fitted, which = 1:npar, alphamax = 0.01, maxpts = 100, 
            delta.t = cutoff/5) 
{
  f.summary <- fitted$summary
  if (is.null(f.summary)) 
    f.summary <- summary(fitted)
  f.model <- fitted$model
  if (is.null(f.model)) 
    stop("Must include model with glm object (run glm with model=T)")
  X <- model.matrix(fitted$terms, f.model, NULL)
  Y <- model.extract(f.model, response)
  f.family <- eval(fitted$call$family)
  std.err <- f.summary$coefficients[, "Std. Error"]
  npar <- length(std.err)
                                        #	class(fitted) <- NULL
  f.class <- class(fitted)
  formula <- fitted$formula
  asgn <- fitted$assign
  p <- .pars <- fitted$coefficients
  pn <- names(p)
  npar <- ncol(X)
  nobs <- nrow(X)
  if (any(which > npar) || any(which < 1)) 
    stop(paste("which must be in the range 1:", npar))
  S.hat <- fitted$deviance
  cutoff <- sqrt(npar * qf(1 - alphamax, npar, nobs - npar))
  s.hat <- sqrt(f.summary$dispersion)
  outmat <- array(0, c(nobs, npar))
  out <- list()
  for (par in which) {
    sgn <- -1
    count <- 1
    varying <- rep(T, npar)
    varying[par] <- F
    .pars <- p
    tau <- double(2 * maxpts)
    par.vals <- array(0, c(2 * maxpts, npar), list(NULL, 
                                                   pn))
    tau[1] <- 0
    par.vals[1, ] <- .pars
    base <- .pars[par]
    profile.par.inc <- delta.t * std.err[par]
    .pars[par] <- .pars[par] - profile.par.inc
    x <- X[, -par, drop = F]
    while (count <= maxpts) {
      if (abs(.pars[par] - base)/std.err[par] > 10 * cutoff) 
        break
      count <- count + 1
      z <- glm.fit(x, Y, family = f.family, offset = .pars[par] * 
                   X[, par])
      tau[count] <- (sgn * sqrt(z$deviance - S.hat))/s.hat
      .pars[-par] <- z$coefficients
      par.vals[count, ] <- .pars
      if (abs(tau[count]) > cutoff) 
        break
      .pars <- .pars + ((.pars - par.vals[count - 1, ]) * 
                        delta.t)/abs(tau[count] - tau[count - 1])
    }
    tau[1:count] <- tau[count:1]
    par.vals[1:count, ] <- par.vals[count:1, ]
    sgn <- 1
    newmax <- count + maxpts
    .pars <- par.vals[count, ]
    while (count <= newmax) {
      .pars <- .pars + ((.pars - par.vals[count - 1, ]) * 
                        delta.t)/abs(tau[count] - tau[count - 1])
      if (abs(.pars[par] - base)/std.err[par] > 10 * cutoff) 
        break
      count <- count + 1
      z <- glm.fit(x, Y, family = f.family, offset = .pars[par] * 
                   X[, par])
      tau[count] <- (sgn * sqrt(z$deviance - S.hat))/s.hat
      .pars[-par] <- z$coefficients
      par.vals[count, ] <- .pars
      if (abs(tau[count]) > cutoff) 
        break
    }
    out[[par]] <- structure(list(tau = tau[1:count], par.vals = par.vals[1:count, 
                                                       ]), class = "data.frame", row.names = as.character(1:count), 
                            parameters = list(par = par, std.err = std.err[par]))
  }
  names(out)[which] <- pn[which]
  orig <- fitted
  class(orig) <- f.class
  attr(out, "original.fit") <- orig
  attr(out, "summary") <- f.summary
  class(out) <- c("profile.glm", "profile.nls", "profile")
  out
}
require(nls)
