.packageName <- "MNP"
mnp <- function(formula, data = parent.frame(), choiceX = NULL,
                cXnames = NULL, base = NULL, latent = FALSE,
                n.draws = 5000, p.var = "Inf", p.df = n.dim+1,
                p.scale = 1, coef.start = 0, cov.start = 1,
                burnin = 0, thin = 0, verbose = FALSE) {   
  call <- match.call()
  mf <- match.call(expand = FALSE)
  mf$choiceX <- mf$cXnames <- mf$base <- mf$n.draws <- mf$latent <-
    mf$p.var <- mf$p.df <- mf$p.scale <- mf$coef.start <-
      mf$cov.start <- mf$verbose <- mf$burnin <- mf$thin <- NULL   
  mf[[1]] <- as.name("model.frame")
  mf$na.action <- 'na.pass'
  mf <- eval.parent(mf)

  ## fix this parameter
  p.alpha0 <- 1
  
  ## obtaining Y
  tmp <- ymatrix.mnp(mf, base=base, extra=TRUE, verbose=verbose)
  Y <- tmp$Y
  MoP <- tmp$MoP
  lev <- tmp$lev
  base <- tmp$base
  p <- tmp$p
  n.dim <- p - 1
  if(verbose)
    cat("\nThe base category is `", base, "'.\n\n", sep="") 
  if (p < 3)
    stop("The number of alternatives should be at least 3.")
  if(verbose) 
    cat("The total number of alternatives is ", p, ".\n\n", sep="") 
  
  ### obtaining X
  tmp <- xmatrix.mnp(formula, data=eval.parent(data),
                     choiceX=call$choiceX, cXnames=cXnames, 
                     base=base, n.dim=n.dim, lev=lev, MoP=MoP,
                     verbose=verbose, extra=TRUE)
  X <- tmp$X
  coefnames <- tmp$coefnames
  n.cov <- ncol(X) / n.dim
  n.obs <- nrow(X)
  if (verbose)
    cat("The dimension of beta is ", n.cov, ".\n\n", sep="")

  ## checking the prior for beta
  p.imp <- FALSE 
  if (p.var == Inf) {
    p.imp <- TRUE
    p.prec <- diag(0, n.cov)
    if (verbose)
      cat("Improper prior will be used for beta.\n\n")
  }
  else if (is.matrix(p.var)) {
    if (ncol(p.var) != n.cov || nrow(p.var) != n.cov)
      stop("The dimension of `p.var' should be ", n.cov, " x ", n.cov, sep="")
    if (sum(sign(eigen(p.var)$values) < 1) > 0)
      stop("`p.var' must be positive definite.")
    p.prec <- solve(p.var)
  }
  else {
    p.var <- diag(p.var, n.cov)
    p.prec <- solve(p.var)
  }
  p.mean <- rep(0, n.cov)

  ## checking prior for Sigma
  p.df <- eval(p.df)
  if (length(p.df) > 1)
    stop("`p.df' must be a positive integer.")
  if (p.df < n.dim)
    stop(paste("`p.df' must be at least ", n.dim, ".", sep=""))
  if (abs(as.integer(p.df) - p.df) > 0)
    stop("`p.df' must be a positive integer.")
  if (!is.matrix(p.scale))  
    p.scale <- diag(p.scale, n.dim)
  if (ncol(p.scale) != n.dim || nrow(p.scale) != n.dim)
    stop("`p.scale' must be ", n.dim, " x ", n.dim, sep="")
  if (sum(sign(eigen(p.scale)$values) < 1) > 0)
    stop("`p.scale' must be positive definite.")
  else if (p.scale[1,1] != 1) {
    p.scale[1,1] <- 1
    warning("p.scale[1,1] will be set to 1.")
  }
  Signames <- NULL
  for(j in 1:n.dim)
    for(k in 1:n.dim)
      if (j<=k)
        Signames <- c(Signames, paste(if(MoP) lev[j] else lev[j+1],
                                      ":", if(MoP) lev[k] else lev[k+1], sep="")) 

  ## checking starting values
  if (length(coef.start) == 1)
    coef.start <- rep(coef.start, n.cov)
  else if (length(coef.start) != n.cov)
    stop(paste("The dimenstion of `coef.start' must be  ",
               n.cov, ".", sep=""))
  if (!is.matrix(cov.start)) {
    cov.start <- diag(n.dim)*cov.start
    cov.start[1,1] <- 1
  }
  else if (ncol(cov.start) != n.dim || nrow(cov.start) != n.dim)
    stop("The dimension of `cov.start' must be ", n.dim, " x ", n.dim, sep="")
  else if (sum(sign(eigen(cov.start)$values) < 1) > 0)
    stop("`cov.start' must be a positive definite matrix.")
  else if (cov.start[1,1] != 1) {
    cov.start[1,1] <- 1
    warning("cov.start[1,1] will be set to 1.")
  }
  
  ## checking thinnig and burnin intervals
  if (burnin < 0)
    stop("`burnin' should be a non-negative integer.") 
  if (thin < 0)
    stop("`thin' should be a non-negative integer.")
  keep <- thin + 1
  
  ## running the algorithm
  if (latent)
    n.par <- n.cov + n.dim*(n.dim+1)/2 + n.dim*n.obs
  else
    n.par <- n.cov + n.dim*(n.dim+1)/2
  if(verbose)
    cat("Starting Gibbs sampler...\n")
  # recoding NA into -1
  Y[is.na(Y)] <- -1

  param <- .C("cMNPgibbs", as.integer(n.dim),
              as.integer(n.cov), as.integer(n.obs), as.integer(n.draws),
              as.double(p.mean), as.double(p.prec), as.integer(p.df),
              as.double(p.scale*p.alpha0), as.double(X), as.integer(Y), 
              as.double(coef.start), as.double(cov.start), 
              as.integer(p.imp), as.integer(burnin), as.integer(keep), 
              as.integer(verbose), as.integer(MoP), as.integer(latent),
              pdStore = double(n.par*floor((n.draws-burnin)/keep)),
              PACKAGE="MNP")$pdStore 
  param <- matrix(param, ncol = n.par,
                  nrow = floor((n.draws-burnin)/keep), byrow=TRUE)
  if (latent) {
    W <- array(as.vector(t(param[,(n.par-n.dim*n.obs+1):n.par])),
               dim = c(n.dim, n.obs, floor((n.draws-burnin)/keep)),
               dimnames = list(lev[-1], rownames(Y), NULL))
    param <- param[,1:(n.par-n.dim*n.obs)]
    }
  else
    W <- NULL
  colnames(param) <- c(coefnames, Signames)
    
  ##recoding -1 back into NA
  Y[Y==-1] <- NA

  ## returning the object
  res <- list(param = param, x = X, y = Y, w = W, call = call, alt = lev,
              n.alt = p, base = base, p.mean = if(p.imp) NULL else p.mean,
              p.var = p.var,
              p.df = p.df, p.scale = p.scale, burnin = burnin, thin = thin) 
  class(res) <- "mnp"
  return(res)
}
  


".onLoad" <- function(lib, pkg) {
  mylib <- dirname(system.file(package = "MNP"))
  title <- paste("MNP:", packageDescription("MNP", lib = mylib)$Title)
  ver <- packageDescription("MNP", lib = mylib)$Version
  url <- packageDescription("MNP", lib = mylib)$URL
  cat(title, "\nVersion:", ver, "\nURL:", url, "\n")
}

predict.mnp <- function(object, newdata = NULL, newdraw = NULL,
                        type = c("prob", "choice", "order", "latent"),
                        verbose = FALSE, ...){

  if (NA %in% match(type, c("prob", "choice", "order", "latent")))
    stop("Invalid input for `type'.")
      
  p <- object$n.alt
  if (is.null(newdraw))
    param <- object$param
  else
    param <- newdraw
  n.cov <- ncol(param) - p*(p-1)/2
  n.draws <- nrow(param)
  coef <- param[,1:n.cov]
  cov <- param[,(n.cov+1):ncol(param)]
  
  ## get X matrix ready
  if (is.null(newdata)) 
    x <- object$x
  else {
    call <- object$call
    x <- xmatrix.mnp(as.formula(call$formula), data = newdata,
                     choiceX = call$choiceX,
                     cXnames = eval(call$cXnames),
                     base = object$base, n.dim = p-1,
                     lev = object$alt, MoP = is.matrix(object$y),
                     verbose = FALSE, extra = FALSE)    
  }

  n.obs <- nrow(x)
  if (verbose)
    cat("There are", n.obs, "observations to predict. Please wait...\n")

  alt <- object$alt
  if (object$base != alt[1]) 
    alt <- c(object$base, alt[1:(length(alt)-1)])
  
  ## computing W
  W <- array(NA, dim=c(p-1, n.obs, n.draws), dimnames=c(alt[2:p],
                                               NULL, 1:n.draws))
  tmp <- floor(n.draws/10)
  inc <- 1
  for (i in 1:n.draws) {
    Sigma <- matrix(0, p-1, p-1)
    count <- 1
    for (j in 1:(p-1)) {
      Sigma[j,j:(p-1)] <- cov[i,count:(count+p-j-1)]
      count <- count + p - j
    }
    diag(Sigma) <- diag(Sigma)/2
    Sigma <- Sigma + t(Sigma)
    for (j in 1:n.obs) 
      W[,j,i] <- matrix(x[j,], ncol=n.cov) %*% matrix(coef[i,]) +
        mvrnorm(1, mu = rep(0, p-1), Sigma = Sigma)
    if (i == inc*tmp & verbose) {
      cat("", inc*10, "percent done.\n")
      inc <- inc + 1
    }
  }
  ans <- list()
  if ("latent" %in% type)
    ans$w <- W
  else
    ans$w <- NULL

  ## computing Y
  Y <- matrix(NA, nrow = n.obs, ncol = n.draws, dimnames=list(NULL, 1:n.draws))
  O <- array(NA, dim=c(p, n.obs, n.draws), dimnames=list(alt, NULL, 1:n.draws))
  for (i in 1:n.obs) 
    for (j in 1:n.draws) {
      Y[i,j] <- alt[match(max(c(0, W[,i,j])), c(0,W[,i,j]))]
      O[,i,j] <- rank(c(0, -W[,i,j]))
    }
  if ("choice" %in% type)
    ans$y <- Y
  else
    ans$y <- NULL
  if ("order" %in% type)
    ans$o <- O
  else
    ans$o <- NULL

  ## computing P
  if ("prob" %in% type) {
    P <- matrix(NA, nrow = n.obs, ncol = p)
    colnames(P) <- alt
    for (i in 1:p)
      P[,i] <- apply(Y==alt[i], 1, mean) 
    ans$p <- P
  }
  else
    ans$p <- NULL

  return(ans)
}
print.mnp <- function (x, digits = max(3, getOption("digits") - 3), ...)
  {
    cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
    param <- apply(x$param, 2, mean)
    if (length(param)) {
      cat("Parameters:\n")
      print.default(format(param, digits = digits), print.gap = 2,
                    quote = FALSE)
    }
    else cat("No parameters\n")
    cat("\n")
    invisible(x)
  }
print.summary.mnp <- function(x, digits = max(3, getOption("digits") - 3), ...) {

  cat("\nCall:\n")
  cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n\n",
      sep = "") 

  cat("\nCoefficients:\n")
  printCoefmat(x$coef.table, digits = digits, na.print = "NA", ...)
  
  cat("\nCovariances:\n")
  printCoefmat(x$cov.table, digits = digits, na.print = "NA", ...)
  
  cat("\nBase category:", x$base)  
  cat("\nNumber of alternatives:", x$n.alt)
  cat("\nNumber of observations:", x$n.obs)
  cat("\nNumber of stored MCMC draws:", x$n.draws)
  cat("\n\n")
  invisible(x)
}
summary.mnp <- function(object, CI=c(2.5, 97.5),...){

  p <- object$n.alt
  param <- object$param
  n.cov <- ncol(param) - p*(p-1)/2
  n.draws <- nrow(param)
  param.table <- cbind(apply(param, 2, mean), apply(param, 2, sd),
                       apply(param, 2, quantile, min(CI)/100),
                       apply(param, 2, quantile, max(CI)/100)) 
  colnames(param.table) <- c("mean", "std.dev.", paste(min(CI), "%", sep=""),
                             paste(max(CI), "%", sep=""))
  
  ans <- list(call=object$call, base = object$base, n.alt=p, n.obs=if(is.matrix(object$y))
              nrow(object$y) else length(object$y), n.draws = n.draws,
              coef.table=param.table[1:n.cov,],
              cov.table=param.table[(n.cov+1):ncol(param),])  
  class(ans) <- "summary.mnp"
  return(ans)
}
xmatrix.mnp <- function(formula, data = parent.frame(), choiceX=NULL,
                        cXnames=NULL, base=NULL, n.dim, lev,
                        MoP=FALSE, verbose=FALSE, extra=FALSE) {
  call <- match.call()
  mf <- match.call(expand = FALSE)
  mf$choiceX <- mf$cXnames <- mf$base <- mf$n.dim <- mf$lev <-
    mf$MoP <- mf$verbose <- mf$extra <- NULL  
  
  ## get variables
  mf[[1]] <- as.name("model.frame.default")
  mf$na.action <- 'na.pass'
  mf <- eval.parent(mf)
  Terms <- attr(mf, "terms")
  X <- model.matrix.default(Terms, mf)
  xvars <- as.character(attr(Terms, "variables"))[-1]
  if ((yvar <- attr(Terms, "response")) > 0)
    xvars <- xvars[-yvar]
  xlev <- if (length(xvars) > 0) {
    xlev <- lapply(mf[xvars], levels)
    xlev[!sapply(xlev, is.null)]
  }
  p <- n.dim + 1
  n.obs <- nrow(X)
  n.cov <- ncol(X)
  
  ## expanding X
  Xcnames <- colnames(X)
  allvnames <- Xnew <- NULL
  for (i in 1:n.cov) {
    Xv <- X[, Xcnames[i]]
    Xtmp <- varnames <- NULL
    for (j in 1:n.dim) {
      allvnames <- c(allvnames, paste(Xcnames[i], ":", if(MoP)
                                      lev[j] else lev[j+1], sep=""))
      for (k in 1:n.dim)
        varnames <- c(varnames, paste(Xcnames[i], ":", if(MoP) lev[j]
        else lev[j+1], sep=""))
      tmp <- matrix(0, nrow = n.obs, ncol = n.dim)
      tmp[, j] <- Xv
      Xtmp <- cbind(Xtmp, tmp)
    }
    colnames(Xtmp) <- varnames
    Xnew <- cbind(Xnew, Xtmp)
  }

  ## checking and adding choice-specific variables
  if (!is.null(choiceX)) {
    cX <- eval(choiceX, data)
    cXn <- unique(names(cX))
    if (sum(is.na(pmatch(cXn, lev))) > 0)
      stop(paste("Error: Invalid input for `choiceX.'\n Some variables do not exist."))
    if(MoP) 
      xbase <- as.matrix(cX[[lev[p]]])
    else if (is.na(pmatch(base, cXn)))
      xbase <- NULL
    else
      xbase <- as.matrix(cX[[base]])
    if (length(cXn) < n.dim)
      stop(paste("Error: Invalid input for `choiceX.'\n You must specify the choice-specific varaibles at least for all non-base categories."))
    if (!is.null(xbase) && length(cXn) != p)
      stop(paste("Error: Invalid input for `choiceX.'\n You must specify the choice-specific variables at least for all non-base categories."))
    if(!is.null(xbase) && verbose)
      cat("The choice-specific variables of the base category are subtracted from the corresponding variables of the non-base categories.\n\n")
    for (i in 1:length(cXnames)) 
      for (j in 1:n.dim) {
        if (length(cXnames) != ncol(as.matrix(cX[[if(MoP) lev[j] else lev[j+1]]])))
            stop(paste("Error: The number of variables in `choiceX' and `cXnames' does not match."))  
        tmp <- matrix(as.matrix(cX[[if(MoP) lev[j] else lev[j+1]]])[,i], ncol=1)
        if (!is.null(xbase)) 
          tmp <- tmp - xbase[,i]
        colnames(tmp) <- paste(cXnames[i], ":", if(MoP) lev[j] else lev[j+1], sep="") 
        Xnew <- cbind(Xnew, tmp)
      }
  }
  if(extra)
    return(list(X=Xnew, coefnames=c(allvnames, cXnames)))
  else
    return(Xnew)
}
ymatrix.mnp <- function(data, base=NULL, extra=FALSE, verbose=verbose) { 
  ## checking and formatting Y
  Y <- model.response(data)
  if (is.matrix(Y)) { # Multinomial ordered Probit model
    for (i in 1:nrow(Y))
      Y[i,] <- match(Y[i,], sort(unique(Y[i,]))) - 1
    p <- ncol(Y)
    lev <- colnames(Y)
    MoP <- TRUE
    if(!is.null(base))
      stop("Error: The last column of the response matrix must be the base category.\n No need to specify `base.'") 
    base <- lev[p]
  }
  else { # standard Multinomial Probit model        
    Y <- as.factor(Y)
    lev <- levels(Y)
    if (!is.null(base))
      if (base %in% lev) {
        lev <- c(base, lev[-pmatch(base, lev)])
        levels(Y) <- lev
      }
      else
        stop(paste("Error: `base' does not exist in the response variable.")) 
    base <- lev[1]
    counts <- table(Y)
    if (any(counts == 0)) {
      warning(paste("group(s)", paste(lev[counts == 0], collapse = " "), "are empty"))
      Y <- factor(Y, levels  = lev[counts > 0])
      lev <- lev[counts > 0]
    }
    p <- length(lev)
    Y <- unclass(Y) - 1
    MoP <- FALSE
  }
  if(extra)
    return(list(Y=Y, MoP=MoP, lev=lev, p=p, base=base))
  else
    return(Y)
}
