#############################################################################
##    Copyright (c) 2010-2022 Rune Haubo Bojesen Christensen
##
##    This file is part of the ordinal package for R (*ordinal*)
##
##    *ordinal* is free software: you can redistribute it and/or modify
##    it under the terms of the GNU General Public License as published by
##    the Free Software Foundation, either version 2 of the License, or
##    (at your option) any later version.
##
##    *ordinal* is distributed in the hope that it will be useful,
##    but WITHOUT ANY WARRANTY; without even the implied warranty of
##    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##    GNU General Public License for more details.
##
##    A copy of the GNU General Public License is available at
##    <https://www.r-project.org/Licenses/> and/or
##    <http://www.gnu.org/licenses/>.
#############################################################################
## This file contains:
## Implementation of various methods for clmm objects.

formatVC <- function(varc, digits = max(3, getOption("digits") - 2))
### "format()" the 'VarCorr' matrix of the random effects -- for
### show()ing
### Borrowed from lme4/R/lmer.R with minor modifications.
{
    recorr <- lapply(varc, attr, "correlation")
    reStdDev <- lapply(varc, attr, "stddev")
    reLens <- unlist(lapply(reStdDev, length))
    nr <- sum(reLens)
    reMat <- array('', c(nr, 4),
		   list(rep.int('', nr),
			c("Groups", "Name", "Variance", "Std.Dev.")))
    reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
    reMat[,2] <- unlist(lapply(varc, colnames))
    reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
    reMat[,4] <- format(unlist(reStdDev), digits = digits)
    if(any(reLens > 1)) {
	maxlen <- max(reLens)
	corr <- do.call("rbind", lapply(recorr, function(x, maxlen)
                                    {
                                        if(is.null(x)) return("")
                                        x <- as(x, "matrix")
                                        cc <- format(round(x, 3), nsmall = 3)
                                        cc[!lower.tri(cc)] <- ""
                                        nr <- dim(cc)[1]
                                        if (nr >= maxlen) return(cc)
                                        cbind(cc, matrix("", nr, maxlen-nr))
                                    },
                                        maxlen))
	colnames(corr) <- c("Corr", rep.int("", maxlen - 1))
        cbind(reMat, corr)
    } else reMat
}


varcov <-
    function(object, format=FALSE,
             digits=max(3, getOption("digits") - 2), ...)
### VarCorr method for model environments - should be the same for
### fitted model objects.
{
    ## Compute variance-covariance matrices of the random effects.
    res <- lapply(object$ST, function(st) {
        ## Variance-covariance matrix for the random effects:
        VC <- tcrossprod(st)
        ## Standard deviations:
        stddev <- sqrt(diag(VC))
        corr <- t(VC / stddev)/stddev
        attr(VC, "stddev") <- stddev
        ## correlation:
        if(NCOL(st) > 1) {
            diag(corr) <- 1
            attr(VC, "correlation") <- corr
        }
        VC
    })
    names(res) <- names(object$dims$nlev.re)
    if(format) noquote(formatVC(res, digits=digits)) else res
}

# VarCorr <- function(x, ...) UseMethod("VarCorr")
VarCorr.clmm <- function(x, ...) varcov(x, ...)

print.clmm <-
  function(x, digits = max(3, getOption("digits") - 3), ...)
{
  if(x$nAGQ >= 2)
    cat(paste("Cumulative Link Mixed Model fitted with the adaptive",
              "Gauss-Hermite \nquadrature approximation with",
              x$nAGQ ,"quadrature points"), "\n\n")
  else if(x$nAGQ <= -1)
    cat(paste("Cumulative Link Mixed Model fitted with the",
              "Gauss-Hermite \nquadrature approximation with",
              abs(x$nAGQ) ,"quadrature points"), "\n\n")
  else
    cat("Cumulative Link Mixed Model fitted with the Laplace approximation\n",
      fill=TRUE)
  cat("formula:", deparse(x$formula), fill=TRUE)
  if(!is.null(data.name <- x$call$data))
    cat("data:   ", deparse(data.name), fill=TRUE)
  if(!is.null(x$call$subset))
    cat("subset: ", deparse(x$call$subset), fill=TRUE)
  cat("\n")

  print(x$info, row.names=FALSE, right=FALSE)

  cat("\nRandom effects:\n")
  print(formatVC(varcov(x), digits=digits), quote=FALSE, ...)
  nlev.char <- paste(names(x$dims$nlev.gf), " ", x$dims$nlev.gf, sep="", collapse=",  ")
  cat("Number of groups: ", nlev.char, "\n")

  if(length(x$beta)) {
    cat("\nCoefficients:\n")
    print(x$beta, digits=digits, ...)
  } else {
    cat("\nNo Coefficients\n")
  }
  if(length(x$alpha) > 0) {
    cat("\nThresholds:\n")
    print(x$alpha, digits=digits, ...)
  }

  if(nzchar(mess <- naprint(x$na.action))) cat("(", mess, ")\n", sep="")
  return(invisible(x))
}

vcov.clmm <- function(object, ...) vcov.clm(object, method="Cholesky")

summary.clmm <- function(object, correlation = FALSE, ...)
{
  if(is.null(object$Hessian))
    stop("Model needs to be fitted with Hess = TRUE")

  nfepar <- object$dims$nfepar
  coef <- matrix(0, nfepar, 4,
                 dimnames = list(names(object$coefficients[1:nfepar]),
                   c("Estimate", "Std. Error", "z value", "Pr(>|z|)")))
  coef[, 1] <- object$coefficients[1:nfepar]
  vc <- try(vcov(object), silent = TRUE)
  if(inherits(vc, "try-error")) {
    warning("Variance-covariance matrix of the parameters is not defined")
    coef[, 2:4] <- NaN
    if(correlation) warning("Correlation matrix is unavailable")
    object$condHess <- NaN
  }
  else {
    coef[, 2] <- sd <- sqrt(diag(vc)[1:nfepar])
    ## Cond is Inf if Hessian contains NaNs:
    object$condHess <-
      if(any(is.na(object$Hessian))) Inf
      else with(eigen(object$Hessian, only.values = TRUE),
                abs(max(values) / min(values)))
    coef[, 3] <- coef[, 1]/coef[, 2]
    coef[, 4] <- 2 * pnorm(abs(coef[, 3]), lower.tail=FALSE)
    if(correlation) ## {
      ## sd <- sqrt(diag(vc))
      object$correlation <- cov2cor(vc)
    ## (vc / sd) / rep(sd, rep(object$edf, object$edf))
  }
  object$info$cond.H <- formatC(object$condHess, digits=1, format="e")
  object$coefficients <- coef
  class(object) <- "summary.clmm"
  return(object)
}

print.summary.clmm <-
  function(x, digits = max(3, getOption("digits") - 3),
           signif.stars = getOption("show.signif.stars"), ...)
{
  if(x$nAGQ >= 2)
    cat(paste("Cumulative Link Mixed Model fitted with the adaptive",
              "Gauss-Hermite \nquadrature approximation with",
              x$nAGQ ,"quadrature points"), "\n\n")
  else if(x$nAGQ <= -1)
    cat(paste("Cumulative Link Mixed Model fitted with the",
              "Gauss-Hermite \nquadrature approximation with",
              abs(x$nAGQ) ,"quadrature points"), "\n\n")
  else
    cat("Cumulative Link Mixed Model fitted with the Laplace approximation\n",
      fill=TRUE)
  cat("formula:", deparse(x$formula), fill=TRUE)
  if(!is.null(data.name <- x$call$data))
    cat("data:   ", deparse(data.name), fill=TRUE)
  if(!is.null(x$call$subset))
    cat("subset: ", deparse(x$call$subset), fill=TRUE)
  cat("\n")

  print(x$info, row.names=FALSE, right=FALSE)

  cat("\nRandom effects:\n")
  print(formatVC(varcov(x), digits=digits), quote=FALSE, ...)
  nlev.char <- paste(names(x$dims$nlev.gf), " ", x$dims$nlev.gf, sep="", collapse=",  ")
  cat("Number of groups: ", nlev.char, "\n")

  nbeta <- length(x$beta)
  nalpha <- length(x$alpha)
  if(nbeta > 0) {
    cat("\nCoefficients:\n")
    printCoefmat(x$coefficients[nalpha + 1:nbeta, , drop=FALSE],
                 digits=digits, signif.stars=signif.stars,
                 has.Pvalue=TRUE, ...)
  } else {
    cat("\nNo Coefficients\n")
  }
  if(nalpha > 0) { ## always true
    cat("\nThreshold coefficients:\n")
    printCoefmat(x$coefficients[seq_len(nalpha), -4, drop=FALSE],
                 digits=digits, has.Pvalue=FALSE, signif.stars=FALSE,
                 ...)
  }

  if(nzchar(mess <- naprint(x$na.action))) cat("(", mess, ")\n", sep="")
  if(!is.null(correl <- x$correlation)) {
    cat("\nCorrelation of Coefficients:\n")
    ll <- lower.tri(correl)
    correl[ll] <- format(round(correl[ll], digits))
    correl[!ll] <- ""
    print(correl[-1, -ncol(correl)], quote = FALSE, ...)
  }
  return(invisible(x))
}

logLik.clmm <- function(object, ...)
  structure(object$logLik, df = object$edf, class = "logLik")

extractAIC.clmm <- function(fit, scale = 0, k = 2, ...) {
  edf <- fit$edf
  c(edf, -2*fit$logLik + k * edf)
}

nobs.clmm <- function(object, ...) object$dims$nobs

edf.clmm <- function(object, ...) object$dims$edf

## anova.clmm <- function(object, ...)
##   anova.clm(object, ...)

anova.clmm <- function(object, ...) {
### This essentially calls anova.clm(object, ...), but the names of
### the models were not displayed correctly in the printed output
### unless the following dodge is enforced.
  mc <- match.call()
  args <- as.list(mc)
  Call <- as.call(c(list(quote(anova.clm)), args[-1]))
  ff <- environment(formula(object))
  pf <- parent.frame()  ## save parent frame in case we need it
  sf <- sys.frames()[[1]]
  ff2 <- environment(object)
  res <- tryCatch(eval(Call, envir=pf),
                  error=function(e) {
                    tryCatch(eval(Call, envir=ff),
                             error=function(e) {
                               tryCatch(eval(Call, envir=ff2),
                                        error=function(e) {
                                          tryCatch(eval(Call, envir=sf),
                                                   error=function(e) {
                                                     "error" })})})})
  if((is.character(res) && res == "error"))
    stop("Unable to evaluate models.")
  res
}


logLik.clmm <- function(object, ...)
  structure(object$logLik, df = object$edf, class = "logLik")

extractAIC.clmm <- function(fit, scale = 0, k = 2, ...) {
  edf <- fit$edf
  c(edf, -2*fit$logLik + k * edf)
}

model.matrix.clmm <- function(object, type = c("design", "B"), ...) {
  type <- match.arg(type)
  mf <- try(model.frame(object), silent=TRUE)
  if(inherits(mf, "try-error"))
    stop("Cannot extract model.matrix: refit model with 'model=TRUE'?")
  if(type == "design") {
    Terms <- terms(object)
    ans <- model.matrix(Terms, mf)
  } else { ## if type == "B":
    stop("type = 'B' not yet implemented")
  }
  return(ans)
}