R/corStruct.R

In nlme: Linear and Nonlinear Mixed Effects Models

###              Classes of correlation structures
###
### Copyright 2005-2020  The R Core team
### Copyright 1997-2003  Jose C. Pinheiro,
###                      Douglas M. Bates <bates@stat.wisc.edu>

#
#  This program 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.
#
#  This program 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
#  http://www.r-project.org/Licenses/
#

### FIXME: For larger problems, the approach of corFactor() / corMatrix() of
### ------ of size sumLenSq := sum(len^2), len := table(groups) sucks !
### Conceptually, one could consider using something like  Matrix::sparseVector(.)
### but would have to deal with corresponding C code as well !


##*## Generics that should be implemented for any corStruct class

corFactor <-
  ## extractor for transpose inverse square root factor of corr matrix
  function(object, ...) UseMethod("corFactor")

corMatrix <-
  ## extractor for correlation matrix or the transpose inverse
  ## square root matrix
  function(object, ...) UseMethod("corMatrix")

###*# Constructor
### There is no constructor function for this class (i.e. no function
### called corStruct) because the class is virtual.

## --- We now check 'sumLenSq' when it is first computed!
## .chkLenSq <- function(nn) {
##   if(nn > .Machine$integer.max)
##     stop(gettextf("'sumLenSq' = %g is too large (larger than maximal integer)",
## 		  nn), domain = NA)
##   else
##     nn
## }

###*# Methods for local generics

corFactor.corStruct <- function(object, ...)
{
  if (!is.null(aux <- attr(object, "factor"))) {
    return(aux)
  }
  corD <- Dim(object)
  ## .chkLenSq(corD[["sumLenSq"]])
  val <- .C(corStruct_factList,
	    as.double(unlist(corMatrix(object))),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
            logDet = double(1))[c("factor", "logDet")]
  lD <- val[["logDet"]]
  val <- val[["factor"]]
  attr(val, "logDet") <- lD
  val
}

corMatrix.corStruct <-
  function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
  if (corr) {
    ## Do not know how to calculate the correlation matrix
      stop(gettextf("do not know how to calculate correlation matrix of %s object",
                    dQuote(class(object)[1])), domain = NA)
  } else {
    ## transpose inverse square root
    corD <- Dim(object,
		if(is.list(covariate)) {
		  if (is.null(names(covariate)))
		    names(covariate) <- seq_along(covariate)
		  rep(names(covariate), lengths(covariate))
		} else
		  rep(1, length(covariate)))
    val <- .C(corStruct_factList,
	      as.double(unlist(corMatrix(object, covariate))),
	      as.integer(unlist(corD)),
	      factor = double(corD[["sumLenSq"]]),
	      logDet = double(1))[c("factor", "logDet")]
    lD <- val[["logDet"]]
    val <- val[["factor"]]
    if (corD[["M"]] > 1) {
      val <- split(val, rep(1:corD[["M"]], (corD[["len"]])^2))
      val <- lapply(val, function(el) {
        nel <- round(sqrt(length(el)))
        array(el, c(nel, nel))
      })
      names(val) <- names(corD[["len"]])
      val <- as.list(val)
    } else {
      val <- array(val, c(corD[["N"]], corD[["N"]]))
    }
    attr(val, "logDet") <- lD
    val
  }
}

###*# Methods for standard generics

as.matrix.corStruct <- function(x, ...) corMatrix(x)

coef.corStruct <-
  ## Accessor for constrained or unconstrained parameters of
  ## corStruct objects
  function(object, unconstrained = TRUE, ...)
{
  if (unconstrained) {
    if (is.null(isFix <- attr(object, "fixed"))) {
      stop("\"corStruct\" object must have a \"fixed\" attribute")
    }
    if (isFix) {
      numeric(0)
    } else {
      as.vector(object)
    }
  } else {
      stop(gettextf("do not know how to obtain parameters of %s object",
                    dQuote(class(object)[1])), domain = NA)
  }
}

`coef<-.corStruct` <- function(object, ..., value)
{
  ## Assignment of the unconstrained parameter of corStruct objects
  value <- as.numeric(value)
  if (length(value) != length(object)) {
    stop("cannot change the length of the parameter of a \"corStruct\" object")
  }
  object[] <- value
  ## updating the factor list and logDet, by forcing a recalculation
  attr(object, "factor") <- NULL
  attr(object, "factor") <- corFactor(object)
  attr(object, "logDet") <- NULL
  attr(object, "logDet") <- logDet(object)
  object
}

Dim.corStruct <- function(object, groups, ...)
{
  if (missing(groups)) return(attr(object, "Dim"))
  ugrp <- unique(groups)
  groups <- factor(groups, levels = ugrp)
  len <- table(groups)
  suml2 <- sum(len^2)
  if(suml2 > .Machine$integer.max)
    stop(gettextf(
      "'sumLenSq := sum(table(groups)^2)' = %g is too large.
 Too large or no groups in your correlation structure?",
		  suml2), call. = FALSE, domain = NA)
  list(N = length(groups),
       M = length(len),
       maxLen = max(len),
       sumLenSq = suml2,
       len = len,
       start = match(ugrp, groups) - 1L)
}

formula.corStruct <-
  ## Accessor for the covariate formula
  function(x, ...) eval(attr(x, "formula"))

getCovariate.corStruct <-
  function(object, form = formula(object), data)
{
  if (!missing(form)) {
    form <- formula(object)
    warning("cannot change 'form'")
  }
  if (is.null(covar <- attr(object, "covariate"))) { # need to calculate it
    if (missing(data)) {
      stop("need data to calculate covariate of \"corStruct\" object")
    }
    covForm <- getCovariateFormula(form)
    grps <- if(!is.null(getGroupsFormula(form)))
                getGroups(object, data = data) ## else NULL
    if (length(all.vars(covForm)) > 0) { # primary covariate present
      if (is.null(grps)) {
        covar <- getCovariate(data, covForm)
      } else {
        if (all(all.vars(covForm) == sapply(splitFormula(covForm, "+"),
                          function(el) deparse(el[[2]])))) {
          covar <- split(getCovariate(data, covForm), grps)
        } else {
          covar <- lapply(split(data, grps), getCovariate, covForm)
        }
      }
    } else {
      if (is.null(grps)) {
        covar <- 1:nrow(data)
      } else {
	covar <- lapply(split(grps, grps), function(x) seq_along(x))
      }
    }
    if (!is.null(grps)) {
      covar <- as.list(covar)
    }
  }
  covar
}

getGroups.corStruct <-
  function(object, form = formula(object), level, data, sep)
{
  if (is.null(val <- attr(object, "groups"))) { # need to calculate
    if (!missing(data)) {
      if ((grpLev <- length(getGroupsFormula(form, asList = TRUE))) > 0) {
        ## use innermost grouping level
        val <- getGroups(data, form, level = grpLev)
        factor(val, levels = unique(as.character(val)))
      } else {
        rep(1, dim(data)[1])
      }
    } else {
      NULL
    }
  } else {
    val
  }
}

Initialize.corStruct <-
  ## Initializes some attributes of corStruct objects
  function(object, data, ...)
{
  form <- formula(object)
  ## obtaining the groups information, if any
  if (!is.null(getGroupsFormula(form))) {
    attr(object, "groups") <- getGroups(object, form, data = data)
    attr(object, "Dim") <- Dim(object, attr(object, "groups"))
  } else {                              # no groups
    attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data))))
  }
  ## obtaining the covariate(s)
  attr(object, "covariate") <- getCovariate(object, data = data)
  object
}

logDet.corStruct <-
  function(object, covariate = getCovariate(object), ...)
{
  if (!is.null(aux <- attr(object, "logDet"))) {
    return(aux)
  }
  if (is.null(aux <- attr(object, "factor"))) {
    ## getting the transpose sqrt factor
    aux <- corMatrix(object, covariate = covariate, corr = FALSE)
  }
  if (is.null(aux1 <- attr(aux, "logDet"))) {
    ## checking for logDet attribute; if not present, get corr matrix
    aux <- corMatrix(object, covariate)
    sum(log(abs(if(is.list(aux)) unlist(lapply(aux, svd.d)) else svd.d(aux)
                )))/2
  } else {
    -aux1
  }
}

## NB, no "nobs"
logLik.corStruct <- function(object, data, ...) -logDet(object)

needUpdate.corStruct <- function(object) FALSE

print.corStruct <- function(x, ...)
{
  if (length(aux <- coef(x, unconstrained = FALSE)) > 0) {
    cat("Correlation structure of class", class(x)[1], "representing\n")
    print(aux, ...)
  } else {
    cat("Uninitialized correlation structure of class", class(x)[1], "\n")
  }
  invisible(x)
}

print.summary.corStruct <- function(x, ...)
{
  class(x) <- attr(x, "oClass")
  cat(paste0("Correlation Structure: ", attr(x, "structName"), "\n"))
  cat(paste(" Formula:", deparse(formula(x)),"\n"))
  cat(" Parameter estimate(s):\n")
  print(coef(x, unconstrained = FALSE), ...)
  invisible(x)
}


recalc.corStruct <-
  function(object, conLin, ...)
{
  conLin[["Xy"]][] <-
    .C(corStruct_recalc,
       Xy = as.double(conLin[["Xy"]]),
       as.integer(unlist(Dim(object))),
       as.integer(ncol(conLin[["Xy"]])),
       as.double(unlist(corFactor(object))))[["Xy"]]
  conLin[["logLik"]] <- conLin[["logLik"]] + logLik(object)
  conLin
}

summary.corStruct <-
  function(object, structName = class(object)[1], ...)
{
  attr(object, "structName") <- structName
  attr(object, "oClass") <- class(object)
  class(object) <- "summary.corStruct"
  object
}

update.corStruct <-
  function(object, data, ...)
{
  object
}

##*## Classes that substitute for (i.e. inherit from) corStruct

###*# corSymm - general, unstructured correlation

####* Constructor

corSymm <-
  ## Constructor for the corSymm class
  function(value = numeric(0), form = ~ 1, fixed = FALSE)
{
  attr(value, "formula") <- form
  attr(value, "fixed") <- fixed
  class(value) <- c("corSymm", "corStruct")
  value
}

###*# Methods for local generics

corFactor.corSymm <- function(object, ...)
{
  corD <- Dim(object)
  val <- .C(symm_factList,
	    as.double(as.vector(object)),
	    as.integer(unlist(attr(object, "covariate"))),
	    as.integer(attr(object, "maxCov")),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  lD <- val[["logDet"]]
  val <- val[["factor"]]
  attr(val, "logDet") <- lD
  val
}

corMatrix.corSymm <-
  function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
  corD <- Dim(object,
	      if(is.list(covariate)) {
		if (is.null(names(covariate)))
		  names(covariate) <- seq_along(covariate)
		rep(names(covariate), lengths(covariate))
	      } else
		rep(1, length(covariate)))
  if (corr) {
    val <- .C(symm_matList,
	      as.double(as.vector(object)),
	      as.integer(unlist(covariate)),
	      as.integer(attr(object, "maxCov")),
	      as.integer(unlist(corD)),
	      mat = double(corD[["sumLenSq"]]))[["mat"]]
    lD <- NULL
  } else {
    val <- .C(symm_factList,
              as.double(as.vector(object)),
              as.integer(unlist(covariate)),
              as.integer(attr(object, "maxCov")),
              as.integer(unlist(corD)),
              factor = double(corD[["sumLenSq"]]),
              logDet = double(1))[c("factor", "logDet")]
    lD <- val[["logDet"]]
    val <- val[["factor"]]
  }
  if (corD[["M"]] > 1) {
    val <- split(val, rep(1:corD[["M"]], (corD[["len"]])^2))
    val <- lapply(val, function(el) {
      nel <- round(sqrt(length(el)))
      array(el, c(nel, nel))
    })
    names(val) <- names(corD[["len"]])
    val <- as.list(val)
  } else {
    val <- array(val, c(corD[["N"]], corD[["N"]]))
  }
  attr(val, "logDet") <- lD
  val
}

###*# Methods for standard generics

coef.corSymm <- function(object, unconstrained = TRUE, ...)
{
  if (unconstrained) {
    if (attr(object, "fixed")) {
      return(numeric(0))
    } else {
      return(as.vector(object))
    }
  }
  mC <- attr(object, "maxCov")
  .C(symm_fullCorr, as.double(object),
     as.integer(mC), corr = double(round(mC * (mC - 1) / 2)))[["corr"]]
}

`coef<-.corSymm` <- function(object, ..., value)
{
  if (length(value) != length(object)) {
    stop("cannot change the length of the parameter of a \"corSymm\" object")
  }
  object[] <- value
  corD <- attr(object, "Dim")
  ## updating the factor list and logDet
  aux <- .C(symm_factList,
	    as.double(as.vector(object)),
	    as.integer(unlist(getCovariate(object))),
	    as.integer(attr(object, "maxCov")),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  attr(object, "factor") <- aux[["factor"]]
  attr(object, "logDet") <- -aux[["logDet"]]
  object
}

Initialize.corSymm <- function(object, data, ...)
{
  if (!is.null(attr(object, "maxCov"))) {# initialized - nothing to do
    return(object)
  }
  object <- NextMethod()

  covar <- attr(object, "covariate")
  if(!is.list(covar)) covar <- list(covar)
  if (any(unlist(lapply(covar, duplicated)))) {
    stop("covariate must have unique values within groups for \"corSymm\" objects")
  }
  covar <- unlist(covar) - 1
  maxCov <- max(uCov <- unique(covar)) + 1
  if (length(uCov) != maxCov) {
    stop("unique values of the covariate  for \"corSymm\" objects must be a sequence of consecutive integers")
  }
  if (Dim(object)[["M"]] > 1) {
    attr(object, "covariate") <- split(covar, getGroups(object))
  } else {
    attr(object, "covariate") <- covar
  }
  attr(object, "maxCov") <- maxCov
  natPar <- as.vector(object)
  if (length(natPar) > 0) {
    ## parameters assumed in constrained form
    if (length(natPar) != round(maxCov * (maxCov - 1) / 2)) {
      stop("initial value for \"corSymm\" parameters of wrong dimension")
    }
    if (max(abs(natPar)) >= 1) {
      stop("initial values for \"corSymm\" must be between -1 and 1")
    }
    natMat <- diag(maxCov)/2
    natMat[lower.tri(natMat)] <- natPar
    natMat <- t(natMat) + natMat
    ## checking if positive-definite
    if (any(eigen(natMat, symmetric=TRUE, only.values=TRUE)$values <= 0)) {
      stop("initial values for \"corSymm\" do not define a positive-definite correlation structure")
    }
    natMat <- chol(natMat)
    uncPar <- numeric(0)
    for(i in 2:maxCov) {
      aux <- acos(natMat[1:(i-1),i]/sqrt(cumsum(natMat[i:1,i]^2)[i:2]))
      uncPar <- c(uncPar, log(aux/(pi - aux)))
    }
    coef(object) <- uncPar
  } else {				# initializing the parameters
    oldAttr <- attributes(object)
    object <- double(round(maxCov * (maxCov - 1) / 2))
    attributes(object) <- oldAttr
    attr(object, "factor") <- corFactor(object)
    attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
  }
  object
}

print.corSymm <- function(x, ...)
{
  if (length(as.vector(x)) > 0 && !is.null(mC <- attr(x, "maxCov"))) {
    aux <- coef.corSymm(x, unconstrained = FALSE)
    val <- diag(mC)
    dimnames(val) <- list(1:mC, 1:mC)
    val[lower.tri(val)] <- aux
    class(val) <- "correlation"
    cat("Correlation structure of class corSymm representing\n")
    print(val, ...)
  }
  else cat("Uninitialized correlation structure of class corSymm\n")
  invisible(x)
}

print.summary.corSymm <- function(x, ...)
{
  if (length(as.vector(x)) > 0 && !is.null(mC <- attr(x, "maxCov"))) {
    cat("Correlation Structure: General\n")
    cat(paste(" Formula:", deparse(formula(x)),"\n"))
    cat(" Parameter estimate(s):\n")
    val <- diag(mC)
    dimnames(val) <- list(1:mC, 1:mC)
    val[lower.tri(val)] <- coef.corSymm(x, unconstrained = FALSE)
    class(val) <- "correlation"
    print(val, ...)
  } else cat("Uninitialized correlation structure of class corSymm\n")
  invisible(x)
}

recalc.corSymm <- function(object, conLin, ...)
{
  val <-
    .C(symm_recalc,
       Xy = as.double(conLin[["Xy"]]),
       as.integer(unlist(Dim(object))),
       as.integer(ncol(conLin[["Xy"]])),
       as.double(as.vector(object)),
       as.integer(unlist(getCovariate(object))),
       as.integer(attr(object, "maxCov")),
       logLik = double(1))[c("Xy", "logLik")]
  conLin[["Xy"]][] <- val[["Xy"]]
  conLin[["logLik"]] <- conLin[["logLik"]] + val[["logLik"]]
  conLin
}

summary.corSymm <- function(object, structName = "General correlation", ...)
{
  attr(object, "structName") <- structName
  class(object) <- "summary.corSymm"
  object
}

###*# corNatural - general correlation in natural parametrization

####* Constructor

corNatural <- function(value = numeric(0), form = ~ 1, fixed = FALSE)
{
  attr(value, "formula") <- form
  attr(value, "fixed") <- fixed
  class(value) <- c("corNatural", "corStruct")
  value
}

###*# Methods for local generics

corFactor.corNatural <- function(object, ...)
{
  corD <- Dim(object)
  val <- .C(nat_factList,
	    as.double(as.vector(object)),
	    as.integer(unlist(attr(object, "covariate"))),
	    as.integer(attr(object, "maxCov")),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
            logDet = double(1))[c("factor", "logDet")]
  lD <- val[["logDet"]]
  val <- val[["factor"]]
  attr(val, "logDet") <- lD
  val
}

corMatrix.corNatural <-
  function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
  corD <- Dim(object,
	      if(is.list(covariate)) {
		if (is.null(names(covariate)))
		  names(covariate) <- seq_along(covariate)
		rep(names(covariate), lengths(covariate))
	      } else
		rep(1, length(covariate)))
  if (corr) {
    val <- .C(nat_matList,
	      as.double(as.vector(object)),
	      as.integer(unlist(covariate)),
	      as.integer(attr(object, "maxCov")),
	      as.integer(unlist(corD)),
              mat = double(corD[["sumLenSq"]]))[["mat"]]
    lD <- NULL
  } else {
    val <- .C(nat_factList,
              as.double(as.vector(object)),
              as.integer(unlist(covariate)),
              as.integer(attr(object, "maxCov")),
              as.integer(unlist(corD)),
              factor = double(corD[["sumLenSq"]]),
              logDet = double(1))[c("factor", "logDet")]
    lD <- val[["logDet"]]
    val <- val[["factor"]]
  }
  if (corD[["M"]] > 1) {
    val <- split(val, rep(1:corD[["M"]], (corD[["len"]])^2))
    val <- lapply(val, function(el) {
        nel <- round(sqrt(length(el)))
        array(el, c(nel, nel))
    })
    names(val) <- names(corD[["len"]])
    val <- as.list(val)
  } else {
    val <- array(val, c(corD[["N"]], corD[["N"]]))
  }
  attr(val, "logDet") <- lD
  val
}

###*# Methods for standard generics

coef.corNatural <-
  function(object, unconstrained = TRUE, ...)
{
  if (unconstrained) {
    if (attr(object, "fixed")) {
      return(numeric(0))
    } else {
      return(as.vector(object))
    }
  }
  mC <- attr(object, "maxCov")
  val <- .C(nat_fullCorr, as.double(object),
            as.integer(mC), corr = double(round(mC * (mC - 1) / 2)))[["corr"]]
  names(val) <- outer(1:mC, 1:mC, function(x,y) paste0("cor(",y,",",x,")"))[
      lower.tri(diag(mC))]
  val
}

`coef<-.corNatural` <- function(object, ..., value)
{
  if (length(value) != length(object)) {
    stop("cannot change the length of the parameter of a \"corNatural\" object")
  }
  object[] <- value
  corD <- attr(object, "Dim")
  ## updating the factor list and logDet
  aux <- .C(nat_factList,
	    as.double(as.vector(object)),
	    as.integer(unlist(getCovariate(object))),
	    as.integer(attr(object, "maxCov")),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  attr(object, "factor") <- aux[["factor"]]
  attr(object, "logDet") <- -aux[["logDet"]]
  object
}

Initialize.corNatural <- function(object, data, ...)
{
  if (!is.null(attr(object, "maxCov"))) {# initialized - nothing to do
    return(object)
  }
  object <- NextMethod()

  covar <- attr(object, "covariate")
  if(!is.list(covar)) covar <- list(covar)
  if (any(unlist(lapply(covar, duplicated)))) {
    stop("covariate must have unique values within groups for \"corNatural\" objects")
  }
  covar <- unlist(covar) - 1
  maxCov <- max(uCov <- unique(covar)) + 1
  if (length(uCov) != maxCov) {
    stop("unique values of the covariate for \"corNatural\" objects must be a sequence of consecutive integers")
  }
  if (Dim(object)[["M"]] > 1) {
    attr(object, "covariate") <- split(covar, getGroups(object))
  } else {
    attr(object, "covariate") <- covar
  }
  attr(object, "maxCov") <- maxCov
  natPar <- as.vector(object)
  if (length(natPar) > 0) {
    ## parameters assumed in constrained form
    if (length(natPar) != round(maxCov * (maxCov - 1) / 2)) {
      stop("initial value for \"corNatural\" parameters of wrong dimension")
    }
    if (max(abs(natPar)) >= 1) {
      stop("initial values for \"corNatural\" must be between -1 and 1")
    }
    natMat <- diag(maxCov)/2
    natMat[lower.tri(natMat)] <- natPar
    natMat <- t(natMat) + natMat
    ## checking if positive-definite
    if (any(eigen(natMat, symmetric=TRUE, only.values=TRUE)$values <= 0)) {
      stop("initial values for \"corNatural\" do not define a positive-definite correlation structure")
    }
    coef(object) <- log((natPar + 1)/(1 - natPar))
  } else {				# initializing the parameters
    oldAttr <- attributes(object)
    object <- double(round(maxCov * (maxCov - 1) / 2))
    attributes(object) <- oldAttr
    attr(object, "factor") <- corFactor(object)
    attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
  }
  object
}

print.corNatural <- function(x, ...)
{
  if (length(as.vector(x)) > 0 &&
      !is.null(mC <- attr(x, "maxCov"))) {
    aux <- coef(x, FALSE)
    val <- diag(mC)
    dimnames(val) <- list(1:mC, 1:mC)
    val[lower.tri(val)] <- aux
    class(val) <- "correlation"
    cat("Correlation structure of class corNatural representing\n")
    print(val, ...)
  }
  else cat("Uninitialized correlation structure of class corNatural\n")
  invisible(x)
}

print.summary.corNatural <-
  function(x, ...)
{
  if (length(as.vector(x)) > 0 &&
      !is.null(mC <- attr(x, "maxCov"))) {
    cat("Correlation Structure: General\n")
    cat(paste(" Formula:", deparse(formula(x)),"\n"))
    cat(" Parameter estimate(s):\n")
    aux <- coef(x, FALSE)
    val <- diag(mC)
    dimnames(val) <- list(1:mC, 1:mC)
    val[lower.tri(val)] <- aux
    class(val) <- "correlation"
    print(val, ...)
  } else cat("Uninitialized correlation structure of class corNatural\n")
  invisible(x)
}

recalc.corNatural <-
  function(object, conLin, ...)
{
  val <-
    .C(nat_recalc,
       Xy = as.double(conLin[["Xy"]]),
       as.integer(unlist(Dim(object))),
       as.integer(ncol(conLin[["Xy"]])),
       as.double(as.vector(object)),
       as.integer(unlist(getCovariate(object))),
       as.integer(attr(object, "maxCov")),
       logLik = double(1))[c("Xy", "logLik")]
  conLin[["Xy"]][] <- val[["Xy"]]
  conLin[["logLik"]] <- conLin[["logLik"]] + val[["logLik"]]
  conLin
}

summary.corNatural <-
  function(object,
           structName = "General correlation, with natural parametrization",
           ...)
{
  attr(object, "structName") <- structName
  class(object) <- "summary.corNatural"
  object
}

###*# corIdent - independent structure === DEPRECATED

####* Constructor

corIdent <-
  ## Constructor for the corIdent class
  function(form = NULL)
{
  .Deprecated("corAR1(0, *)", "nlme")
  value <- numeric(0)
  attr(value, "formula") <- form
  attr(value, "fixed") <- TRUE
  class(value) <- c("corIdent", "corStruct")
  value
}

###*# Methods for local generics

corMatrix.corIdent <-
  function(object, covariate = getCovariate(object), corr, ...)
{
  .Deprecated(package="nlme", old = "corMatrix.corIdent(): \"corIdent\" class")
  if(is.list(covariate)) {# by group
    as.list(lapply(covariate, function(el, object) corMatrix(object, el)))
  } else {
    diag(length(covariate))
  }
}

###*# Methods for standard generics

coef.corIdent <-
  function(object, unconstrained = TRUE, ...) numeric(0)

`coef<-.corIdent` <- function(object, ..., value) object

Initialize.corIdent <-
  function(object, data, ...)
{
  .Deprecated(package="nlme", old = "Initialize.corIdent(): \"corIdent\" class")
  attr(object, "logDet") <- 0
  object
}

logDet.corIdent <- function(object, covariate, ...) 0

recalc.corIdent <- function(object, conLin, ...) conLin

summary.corIdent <-
  function(object, structName = "Independent", ...)
{
  .Deprecated(package="nlme", old = "summary.corIdent(): \"corIdent\" class")
  summary.corStruct(object, structName)
}

###*# corAR1 - autoregressive of order one structure

####* Constructor

corAR1 <-
  ## Constructor for the corAR1 class
  function(value = 0, form = ~ 1, fixed = FALSE)
{
  if (abs(value) >= 1) {
    stop("parameter in AR(1) structure must be between -1 and 1")
  }
  value <- log((1 + value)/( 1 - value))
  attr(value, "formula") <- form
  attr(value, "fixed") <- fixed
  class(value) <- c("corAR1", "corStruct")
  value
}

###*# Methods for local generics

corFactor.corAR1 <- function(object, ...)
{
  corD <- Dim(object)
  if(corD[["sumLenSq"]] > .Machine$integer.max)
    stop(gettextf("'sumLenSq' = %g is too large (larger than maximal integer)",
		  corD[["sumLenSq"]]), domain = NA)
  val <- .C(AR1_factList,
	    as.double(as.vector(object)),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]), ## of size n^2 -- too large !!
	    logDet = double(1))[c("factor", "logDet")]
  lD <- val[["logDet"]]
  val <- val[["factor"]]
  attr(val, "logDet") <- lD
  val
}

corMatrix.corAR1 <-
  function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
  corD <- Dim(object,
	      if(is.list(covariate)) {
		if (is.null(names(covariate)))
		  names(covariate) <- seq_along(covariate)
		rep(names(covariate), lengths(covariate))
	      } else
		rep(1, length(covariate)))
  if (corr) {
    val <- .C(AR1_matList,
	      as.double(as.vector(object)),
	      as.integer(unlist(corD)),
	      mat = double(corD[["sumLenSq"]]))[["mat"]]
    lD <- NULL
  } else {
    val <- .C(AR1_factList,
              as.double(as.vector(object)),
              as.integer(unlist(corD)),
              factor = double(corD[["sumLenSq"]]),
              logDet = double(1))[c("factor", "logDet")]
    lD <- val[["logDet"]]
    val <- val[["factor"]]
  }
  if (corD[["M"]] > 1) {
    val <- split(val, rep(1:corD[["M"]], (corD[["len"]])^2))
    val <- lapply(val, function(el) {
      nel <- round(sqrt(length(el)))
      array(el, c(nel, nel))
    })
    names(val) <- names(corD[["len"]])
    val <- as.list(val)
  } else {
    val <- array(val, c(corD[["N"]], corD[["N"]]))
  }
  attr(val, "logDet") <- lD
  val
}

###*# Methods for standard generics

coef.corAR1 <-
  function(object, unconstrained = TRUE, ...)
{
  if (unconstrained) {
    if (attr(object, "fixed")) {
      return(numeric(0))
    } else {
      return(as.vector(object))
    }
  }
  aux <- exp(as.vector(object))
  aux <- c((aux - 1)/(aux + 1))
  names(aux) <- "Phi"
  aux
}

`coef<-.corAR1` <- function(object, ..., value)
{
  if (length(value) != length(object)) {
    stop("cannot change the length of the parameter of a \"corAR1\" object")
  }
  object[] <- value
  corD <- attr(object, "Dim")
  ## updating the factor list and logDet
  aux <- .C(AR1_factList,
	    as.double(as.vector(object)),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  attr(object, "factor") <- aux[["factor"]]
  attr(object, "logDet") <- -aux[["logDet"]]
  object
}

Initialize.corAR1 <-
  ## Initializes corAR1 objects
  function(object, data, ...)
{
  object <- NextMethod()
  covar <- attr(object, "covariate")
  if(!is.list(covar)) covar <- list(covar)
  if (any(unlist(lapply(covar, duplicated)))) {
    stop("covariate must have unique values within groups for \"corAR1\" objects")
  }
  if (any(unlist(lapply(covar, diff)) != 1)) {
    ## Cannot use formulas for inverse of square root matrix
    ## will convert to class ARMA(1,0)
    attr(object, "p") <- 1
    attr(object, "q") <- 0
    class(object) <- c("corARMA", "corStruct")
    Initialize(object, data)
  } else {
    ## obtaining the factor list and logDet
    attr(object, "factor") <- corFactor(object)
    attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
    object
  }
}

recalc.corAR1 <-
  function(object, conLin, ...)
{
  val <-
    .C(AR1_recalc,
       Xy = as.double(conLin[["Xy"]]),
       as.integer(unlist(Dim(object))),
       as.integer(ncol(conLin[["Xy"]])),
       as.double(as.vector(object)),
       logLik = double(1))[c("Xy", "logLik")]
  conLin[["Xy"]][] <- val[["Xy"]]
  conLin[["logLik"]] <- conLin[["logLik"]] + val[["logLik"]]
  conLin
}

summary.corAR1 <-
  function(object, structName = "AR(1)", ...)
{
  summary.corStruct(object, structName)
}

####*# corCAR1 - continuous time autoregressive of order one structure

#####* Constructor

corCAR1 <-
  ## Constructor for the corCAR1 class
  function(value = 0.2, form = ~ 1, fixed = FALSE)
{
  if (value <= 0 | value >= 1) {
    stop("parameter in CAR(1) structure must be between 0 and 1")
  }
  value <- log(value / (1 - value))
  attr(value, "formula") <- form
  attr(value, "fixed") <- fixed
  class(value) <- c("corCAR1", "corStruct")
  value
}


###*# Methods for local generics

corFactor.corCAR1 <-
  function(object, ...)
{
  corD <- Dim(object)
  val <- .C(CAR1_factList,
	    as.double(as.vector(object)),
	    as.double(unlist(attr(object, "covariate"))),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  lD <- val[["logDet"]]
  val <- val[["factor"]]
  attr(val, "logDet") <- lD
  val
}

corMatrix.corCAR1 <-
  function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
  corD <- Dim(object,
	      if(is.list(covariate)) {
		if (is.null(names(covariate)))
		  names(covariate) <- seq_along(covariate)
		rep(names(covariate), lengths(covariate))
	      } else
		rep(1, length(covariate)))
  if (corr) {
    val <- .C(CAR1_matList,
	      as.double(as.vector(object)),
	      as.double(unlist(covariate)),
	      as.integer(unlist(corD)),
	      mat = double(corD[["sumLenSq"]]))[["mat"]]
    lD <- NULL
  } else {
    val <- .C(CAR1_factList,
              as.double(as.vector(object)),
              as.double(unlist(covariate)),
              as.integer(unlist(corD)),
              factor = double(corD[["sumLenSq"]]),
              logDet = double(1))[c("factor", "logDet")]
    lD <- val[["logDet"]]
    val <- val[["factor"]]
  }
  if (corD[["M"]] > 1) {
    val <- split(val, rep(1:corD[["M"]], (corD[["len"]])^2))
    val <- lapply(val, function(el) {
      nel <- round(sqrt(length(el)))
      array(el, c(nel, nel))
    })
    names(val) <- names(corD[["len"]])
    val <- as.list(val)
  } else {
    val <- array(val, c(corD[["N"]], corD[["N"]]))
  }
  attr(val, "logDet") <- lD
  val
}

###*# Methods for standard generics

coef.corCAR1 <-
  function(object, unconstrained = TRUE, ...)
{
  if (unconstrained) {
    if (attr(object, "fixed")) {
      return(numeric(0))
    } else {
      return(as.vector(object))
    }
  }
  aux <- c(exp(as.vector(object)))
  aux <- aux/(1+aux)
  names(aux) <- "Phi"
  aux
}

`coef<-.corCAR1` <- function(object, ..., value)
{
  if (length(value) != length(object)) {
    stop("cannot change the length of the parameter of a \"corCAR1\" object")
  }
  object[] <- value
  corD <- attr(object, "Dim")
  ## updating the factor list and logDet
  aux <- .C(CAR1_factList,
	    as.double(as.vector(object)),
	    as.double(unlist(getCovariate(object))),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  attr(object, "factor") <- aux[["factor"]]
  attr(object, "logDet") <- -aux[["logDet"]]
  object
}

Initialize.corCAR1 <-
  ## Initializes corCAR1 objects
  function(object, data, ...)
{
  object <- NextMethod()
  covar <- attr(object, "covariate")
  if(!is.list(covar)) covar <- list(covar)

  if (any(unlist(lapply(covar, duplicated)))) {
    stop("covariate must have unique values within groups for \"corCAR1\" objects")
  }
  attr(object, "factor") <- corFactor(object)
  attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
  object
}

recalc.corCAR1 <-
  function(object, conLin, ...)
{
  val <-
    .C(CAR1_recalc,
     Xy = as.double(conLin[["Xy"]]),
     as.integer(unlist(Dim(object))),
     as.integer(ncol(conLin[["Xy"]])),
     as.double(as.vector(object)),
     as.double(unlist(getCovariate(object))),
     logLik = double(1))[c("Xy", "logLik")]
  conLin[["Xy"]][] <- val[["Xy"]]
  conLin[["logLik"]] <- conLin[["logLik"]] + val[["logLik"]]
  conLin
}

summary.corCAR1 <-
  function(object, structName = "Continuous AR(1)", ...)
{
  summary.corStruct(object, structName)
}

###*# corARMA - autoregressive-moving average structures

####* Constructor

corARMA <-
  ## Constructor for the corARMA class
  function(value = double(p + q), form = ~ 1, p = 0, q = 0, fixed = FALSE)
{
  if (!(p >= 0 && (p == round(p)))) {
    stop("autoregressive order must be a non-negative integer")
  }
  if (!(q >= 0 && (q == round(q)))) {
    stop("moving average order must be a non-negative integer")
  }
  if (0 == (p + q)) {
    return(corIdent())
  }
  if (length(value) != p + q) {
    stop("initial value for parameter of wrong length")
  }
  if (max(abs(value)) >= 1) {
    stop("parameters in ARMA structure must be < 1 in absolute value")
  }
  ## unconstrained parameters
  value <- .C(ARMA_unconstCoef,
	      as.integer(p),
	      as.integer(q),
	      pars = as.double(value))$pars
  attributes(value) <- list(formula = form, p = p, q = q, fixed = fixed)
  class(value) <- c("corARMA", "corStruct")
  value
}


###*# Methods for local generics

corFactor.corARMA <-
  function(object, ...)
{
    maxLag <- attr(object, "maxLag")
    if(is.null(maxLag)) stop("'object' has not been Initialize()d")
  corD <- Dim(object)
  val <- .C(ARMA_factList,
	    as.double(as.vector(object)),
	    as.integer(attr(object, "p")),
	    as.integer(attr(object, "q")),
	    as.integer(unlist(attr(object, "covariate"))),
	    as.integer(attr(object, "maxLag")),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  lD <- val[["logDet"]]
  val <- val[["factor"]]
  attr(val, "logDet") <- lD
  val
}


corMatrix.corARMA <-
  function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
  corD <- Dim(object,
	      if(is.list(covariate)) {
		if (is.null(names(covariate)))
		  names(covariate) <- seq_along(covariate)
		rep(names(covariate), lengths(covariate))
	      } else
		rep(1, length(covariate)))
  p <- attr(object, "p")
  q <- attr(object, "q")
  maxLag <- attr(object, "maxLag")
  if(is.null(maxLag)) stop("'object' has not been Initialize()d")
  if (corr) {
    val <- .C(ARMA_matList,
	      as.double(as.vector(object)),
	      as.integer(p),
	      as.integer(q),
	      as.integer(unlist(covariate)),
	      as.integer(maxLag),
	      as.integer(unlist(corD)),
	      mat = double(corD[["sumLenSq"]]))[["mat"]]
    lD <- NULL
  } else {
    val <- .C(ARMA_factList,
              as.double(as.vector(object)),
              as.integer(attr(object, "p")),
              as.integer(attr(object, "q")),
              as.integer(unlist(covariate)),
              as.integer(attr(object, "maxLag")),
              as.integer(unlist(corD)),
              factor = double(corD[["sumLenSq"]]),
              logDet = double(1))[c("factor", "logDet")]
    lD <- val[["logDet"]]
    val <- val[["factor"]]
  }
  if (corD[["M"]] > 1) {
    val <- split(val, rep(1:corD[["M"]], (corD[["len"]])^2))
    val <- lapply(val, function(el) {
      nel <- round(sqrt(length(el)))
      array(el, c(nel, nel))
    })
    names(val) <- names(corD[["len"]])
    val <- as.list(val)
  } else {
    val <- array(val, c(corD[["N"]], corD[["N"]]))
  }
  attr(val, "logDet") <- lD
  val
}

###*# Methods for standard generics

coef.corARMA <-
  function(object, unconstrained = TRUE, ...)
{
  if (attr(object, "fixed") && unconstrained) {
    return(numeric(0))
  }
  val <-  as.vector(object)
  if (!unconstrained) {
    p <- attr(object, "p")
    q <- attr(object, "q")
    nams <- NULL
    if (p > 0) {
      nams <- paste(rep("Phi", p), 1:p, sep="")
    }
    if (q > 0) {
      nams <- c(nams, paste(rep("Theta", q), 1:q, sep=""))
    }
    val <- c(.C(ARMA_constCoef, as.integer(attr(object,"p")),
		as.integer(attr(object,"q")),
		pars = as.double(val))$pars)
    names(val) <- nams
  }
  val
}

`coef<-.corARMA` <- function(object, ..., value)
{
  maxLag <- attr(object, "maxLag")
  if(is.null(maxLag)) stop("'object' has not been Initialize()d")
  if (length(value) != length(object)) {
    stop("cannot change the length of the parameter of a \"corARMA\" object")
  }
  p <- attr(object, "p")
  q <- attr(object, "q")
  object[] <- value
  ## updating the factor list and logDet
  corD <- Dim(object)
  aux <- .C(ARMA_factList,
	    as.double(as.vector(object)),
	    as.integer(p),
	    as.integer(q),
	    as.integer(unlist(getCovariate(object))),
	    as.integer(attr(object, "maxLag")),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  attr(object, "factor") <- aux[["factor"]]
  attr(object, "logDet") <- -aux[["logDet"]]
  object
}

Initialize.corARMA <-
  function(object, data, ...)
{
  ## Initializes corARMA objects
  object <- NextMethod()
  covar <- attr(object, "covariate")
  if(!is.list(covar)) covar <- list(covar)
  if (any(unlist(lapply(covar, duplicated)))) {
    stop("covariate must have unique values within groups for \"corARMA\" objects")
  }
  if ((attr(object, "p") == 1) && (attr(object, "q") == 0) &&
     all(unlist(lapply(covar, diff)) == 1)) {
    ## Use AR1 methods instead
    class(object) <- c("corAR1", "corStruct")
    Initialize(object, data)
  } else {
    attr(object, "maxLag") <-
      max(unlist(lapply(covar, function(el) max(abs(outer(el,el,"-"))))))
    attr(object, "factor") <- corFactor(object)
    attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
    object
  }
}

recalc.corARMA <-
  function(object, conLin, ...)
{
    maxLag <- attr(object, "maxLag")
    if(is.null(maxLag)) stop("'object' has not been Initialize()d")
  val <-
    .C(ARMA_recalc,
     Xy = as.double(conLin[["Xy"]]),
     as.integer(unlist(Dim(object))),
     as.integer(ncol(conLin[["Xy"]])),
     as.double(as.vector(object)),
     as.integer(attr(object, "p")),
     as.integer(attr(object, "q")),
     as.integer(unlist(getCovariate(object))),
     as.integer(attr(object, "maxLag")),
     logLik = double(1))[c("Xy", "logLik")]
  conLin[["Xy"]][] <- val[["Xy"]]
  conLin[["logLik"]] <- conLin[["logLik"]] + val[["logLik"]]
  conLin
}

summary.corARMA <-
  function(object, structName = paste("ARMA(",attr(object,"p"),",",
		     attr(object,"q"), ")", sep = ""), ...)
{
  summary.corStruct(object, structName)
}

###*# corCompSymm - Compound symmetry structure structure

####* Constructor

corCompSymm <-
  ## Constructor for the corCompSymm class
  function(value = 0, form = ~ 1, fixed = FALSE)
{
  if (abs(value) >= 1) {
    stop("parameter in \"corCompSymm\" structure must be < 1 in absolute value")
  }
  attr(value, "formula") <- form
  attr(value, "fixed") <- fixed
  class(value) <- c("corCompSymm", "corStruct")
  value
}

###*# Methods for local generics

corFactor.corCompSymm <-
  function(object, ...)
{
  corD <- Dim(object)
  val <- .C(compSymm_factList,
	    as.double(as.vector(object)),
	    as.double(attr(object, "inf")),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  lD <- val[["logDet"]]
  val <- val[["factor"]]
  attr(val, "logDet") <- lD
  val
}

corMatrix.corCompSymm <-
  function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
  corD <- Dim(object,
	      if(is.list(covariate)) {
		if (is.null(names(covariate)))
		  names(covariate) <- seq_along(covariate)
		rep(names(covariate), lengths(covariate))
	      } else
		rep(1, length(covariate)))
  if (corr) {
    val <- .C(compSymm_matList,
	      as.double(as.vector(object)),
	      as.double(attr(object, "inf")),
	      as.integer(unlist(corD)),
	      mat = double(corD[["sumLenSq"]]))[["mat"]]
    lD <- NULL
  } else {
    val <- .C(compSymm_factList,
              as.double(as.vector(object)),
              as.double(attr(object, "inf")),
              as.integer(unlist(corD)),
              factor = double(corD[["sumLenSq"]]),
              logDet = double(1))[c("factor", "logDet")]
    lD <- val[["logDet"]]
    val <- val[["factor"]]
  }
  if (corD[["M"]] > 1) {
    val <- split(val, rep(1:corD[["M"]], (corD[["len"]])^2))
    val <- lapply(val, function(el) {
      nel <- round(sqrt(length(el)))
      array(el, c(nel, nel))
    })
    names(val) <- names(corD[["len"]])
    val <- as.list(val)
  } else {
    val <- array(val, c(corD[["N"]], corD[["N"]]))
  }
  attr(val, "logDet") <- lD
  val
}

###*# Methods for local generics

coef.corCompSymm <-
  function(object, unconstrained = TRUE, ...)
{
  if (unconstrained) {
    if (attr(object, "fixed")) {
      return(numeric(0))
    } else {
      return(as.vector(object))
    }
  }
  val <- exp(as.vector(object))
  val <- c((val + attr(object, "inf"))/(val + 1))
  names(val) <- "Rho"
  val
}

`coef<-.corCompSymm` <- function(object, ..., value)
{
  if (length(value) != length(object)) {
    stop("cannot change the length of the parameter of a \"corCompSymm\" object")
  }
  object[] <- value
  corD <- attr(object, "Dim")
  ## updating the factor list and logDet
  aux <- .C(compSymm_factList,
	    as.double(as.vector(object)),
	    as.double(attr(object, "inf")),
	    as.integer(unlist(corD)),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  attr(object, "factor") <- aux[["factor"]]
  attr(object, "logDet") <- -aux[["logDet"]]
  object
}

Initialize.corCompSymm <-
  ## Initializes corCompSymm objects
  function(object, data, ...)
{
  if (!is.null(attr(object, "inf"))) {   # initialized - nothing to do
    return(object)
  }
  object <- NextMethod()
  natPar <- as.vector(object)
  corD <- Dim(object)
  if (natPar <= (attr(object, "inf") <- -1/(corD[["maxLen"]] - 1))) {
      stop(gettextf("initial value in \"corCompSymm\" must be greater than %s",
                    attr(object, "inf")), domain = NA)
  }
  object[] <- log((natPar - attr(object, "inf"))/(1 - natPar))
  attr(object, "factor") <- corFactor(object)
  attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
  object
}

print.corCompSymm <- function(x, ...)
{
  if (length(as.vector(x)) > 0 && !is.null(attr(x, "inf"))) {
    NextMethod()
  } else {
    cat("Uninitialized correlation structure of class corCompSymm\n")
    invisible(x)
  }
}

print.summary.corCompSymm <- print.corCompSymm

recalc.corCompSymm <-
  function(object, conLin, ...)
{
  val <-
    .C(compSymm_recalc,
       Xy = as.double(conLin[["Xy"]]),
       as.integer(unlist(Dim(object))),
       as.integer(ncol(conLin[["Xy"]])),
       as.double(as.vector(object)),
       as.double(attr(object, "inf")),
       logLik = double(1))[c("Xy", "logLik")]
  conLin[["Xy"]][] <- val[["Xy"]]
  conLin[["logLik"]] <- conLin[["logLik"]] + val[["logLik"]]
  conLin
}

summary.corCompSymm <-
  function(object, structName = "Compound symmetry", ...)
{
  object <- summary.corStruct(object, structName)
  class(object) <- c("summary.corCompSymm", class(object))
  object
}

####*# corHF - Huyn-Feldt structure

#corHF <-
#  ## Constructor for the corHuynFeldt class
#  function(value = numeric(0), form = ~ 1)
#{
#  attr(value, "formula") <- form
#  class(value) <- c("corHF", "corStruct")
#  value
#}

####*# Methods for local generics

#corFactor.corHF <-
#  function(object)
#{
#  corD <- Dim(object)
#  val <- .C("HF_factList",
#	    as.double(as.vector(object)),
#	    as.integer(attr(object, "maxCov")),
#	    as.integer(unlist(getCovariate(object))),
#	    as.integer(unlist(corD)),
#	    factor = double(corD[["sumLenSq"]]),
#	    logDet = double(1))[c("factor", "logDet")]
#  lD <- val[["logDet"]]
#  val <- val[["factor"]]
#  attr(val, "logDet") <- lD
#  val
#}

#corMatrix.corHF <-
#  function(object, covariate = getCovariate(object), corr = TRUE)
#{
#   corD <- Dim(object,
# 	      if(is.list(covariate)) {
# 		if (is.null(names(covariate)))
# 		  names(covariate) <- seq_along(covariate)
# 		rep(names(covariate), lengths(covariate))
# 	      } else
# 		rep(1, length(covariate)))
#  if (corr) {
#    val <- .C("HF_matList",
#	      as.double(as.vector(object)),
#	      as.integer(attr(object, "maxCov")),
#	      as.integer(unlist(covariate)),
#	      as.integer(unlist(corD)),
#	      mat = double(corD[["sumLenSq"]]))[["mat"]]
#    lD <- NULL
#  } else {
#    val <- .C("HF_factList",
#              as.double(as.vector(object)),
#              as.integer(attr(object, "maxCov")),
#              as.integer(unlist(covariate)),
#              as.integer(unlist(corD)),
#              factor = double(corD[["sumLenSq"]]),
#              logDet = double(1))[c("factor", "logDet")]
#    lD <- val[["logDet"]]
#    val <- val[["factor"]]
#  }
#  if (corD[["M"]] > 1) {
#    val <- split(val, rep(1:corD[["M"]], (corD[["len"]])^2))
#    val <- lapply(val, function(el) {
#      nel <- round(sqrt(length(el)))
#      array(el, c(nel, nel))
#    })
#    names(val) <- names(corD[["len"]])
#  } else {
#    val <- array(val, c(corD[["N"]], corD[["N"]]))
#  }
#  attr(val, "logDet") <- lD
#  val
#}

####*# Methods for standard generics

#coef.corHF <-
#  function(object, unconstrained = TRUE)
#{
#  aux <- as.vector(object)
#  if (!unconstrained) {
#    aux <- 2 * (exp(aux) + attr(object, "inf")) + 1
#  }
#  aux
#}

#"coef<-.corHF" <-
#  function(object, value)
#{
#  if (length(value) != length(object)) {
#    stop("Cannot change the length of the parameter of a corStruct object")
#  }
#  object[] <- value
#  corD <- attr(object, "Dim")
#  ## updating the factor list and logDet
#  aux <- .C("HF_factList",
#	    as.double(as.vector(object)),
#	    as.integer(attr(object, "maxCov")),
#	    as.integer(unlist(getCovariate(object))),
#	    as.integer(unlist(corD)),
#	    factor = double(corD[["sumLenSq"]]),
#	    logDet = double(1))[c("factor", "logDet")]
#  attr(object, "factor") <- aux[["factor"]]
#  attr(object, "logDet") <- -aux[["logDet"]]
#  object
#}

#initialize.corHF <-
#  function(object, data, ...)
#{
#  if (!is.null(attr(object, "inf"))) {   # initialized - nothing to do
#    return(object)
#  }
#  object <- NextMethod()
#  covar <- attr(object, "covariate")
#  if (is.list(covar)) {
#    attr(object, "covariate") <- covar <-
#      lapply(covar, function(el) el - 1)
#  } else {
#    attr(object, "covariate") <- covar <- covar - 1
#    covar <- list(covar)
#  }
#  if (any(unlist(lapply(covar, duplicated)))) {
#    stop(paste("Covariate must have unique values",
#               "within groups for corHF objects"))
#  }
#  maxCov <- max(uCov <- unique(unlist(covar))) + 1
#  if (length(uCov) != maxCov) {
#    stop(paste("Unique values of the covariate  for \"corHF\"",
#               "objects must be a sequence of consecutive integers"))
#  }
#  attr(object, "maxCov") <- maxCov
#  attr(object, "inf") <- -1/(2*maxCov)
#  natPar <- as.vector(object)
#  if (length(natPar) > 0) {
#    if (length(aux) != attr(object, "maxCov"))
#      stop("Initial value for Huyn-Feldt parameters of wrong dimension")
#    ## verifying if initial values satisfy constraints
#    if (any(natPar <= attr(object, "inf"))) {
#      stop(paste("Initial values for \"corHF\" parameters",
#		 "must be > than", attr(object, "inf")))
#    }
#    object[] <- log(natPar - attr(object, "inf"))
#  } else {				# initializing the parameters
#    oldAttr <- attributes(object)
#    object <- log(rep(-attr(object, "inf"), att(object, "maxCov")))
#    attributes(object) <- oldAttr
#  }
#  attr(object, "factor") <- corFactor(object)
#  attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
#  object
#}

#print.corHF <-
#  function(x, ...)
#{
#  if (length(as.vector(x)) > 0 && !is.null(attr(object, "maxCov")))
#    NextMethod()
#  else cat("Uninitialized correlation structure of class corHF\n")
#}

#recalc.corHF <-
#  function(object, conLin)
#{
#  val <-
#    .C("HF_recalc",
#       Xy = as.double(conLin[["Xy"]]),
#       as.integer(unlist(Dim(object))),
#       as.integer(ncol(conLin[["Xy"]])),
#       as.double(as.vector(object)),
#       as.integer(unlist(getCovariate(object))),
#       as.integer(attr(object, "maxCov")),
#       logLik = double(1))[c("Xy", "logLik")]
#  conLin[["Xy"]][] <- val[["Xy"]]
#  conLin[["logLik"]] <- conLin[["logLik"]] + val[["logLik"]]
#  conLin
#}

#summary.corHF <-
#  function(object, structName = "Huyn-Feldt")
#{
#  summary.corStruct(object, structName)
#}

###*# corSpatial - a virtual class of spatial correlation structures

###*# Constructor

corSpatial <-
  ## Constructor for the corSpatial class
  function(value = numeric(0), form = ~ 1, nugget = FALSE,
	   type = c("spherical", "exponential", "gaussian", "linear",
             "rational"),
	   metric = c("euclidean", "maximum", "manhattan"), fixed = FALSE)
{
  spClass <- c(spherical = "corSpher",
	       exponential = "corExp",
	       gaussian = "corGaus",
	       linear = "corLin",
	       rational = "corRatio")[match.arg(type)]
  structure(value,
	    "formula" = form,
	    "nugget" = nugget,
	    "metric" = match.arg(metric),
	    "fixed" = fixed,
	    class = c(spClass, "corSpatial", "corStruct"))
}

###*# Methods for local generics

corFactor.corSpatial <- function(object, ...)
{
  corD <- Dim(object)
  val <- .C(spatial_factList,
	    as.double(as.vector(object)),
	    as.integer(attr(object, "nugget")),
	    as.double(unlist(getCovariate(object))),
	    as.integer(unlist(corD)),
	    as.double(attr(object, "minD")),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  structure(val[["factor"]], logDet = val[["logDet"]])
}

corMatrix.corSpatial <-
  function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
  nRt <- function(vec) round((1 + sqrt(1 + 8 * length(vec))) / 2)
  corD <- Dim(object,
	      if(is.list(covariate)) {
		if (is.null(names(covariate)))
		  names(covariate) <- seq_along(covariate)
		rep(names(covariate), vapply(covariate, nRt, numeric(1)))
	      }
	      else
		rep(1, nRt(covariate)))
  if (corr) {
    val <- .C(spatial_matList,
	      as.double(as.vector(object)),
	      as.integer(attr(object, "nugget")),
	      as.double(unlist(covariate)),
	      as.integer(unlist(corD)),
	      as.double(attr(object, "minD")),
	      mat = double(corD[["sumLenSq"]]))[["mat"]]
    lD <- NULL
  } else {
    val <- .C(spatial_factList,
              as.double(as.vector(object)),
              as.integer(attr(object, "nugget")),
              as.double(unlist(covariate)),
              as.integer(unlist(corD)),
              as.double(attr(object, "minD")),
              factor = double(corD[["sumLenSq"]]),
              logDet = double(1))[c("factor", "logDet")]
    lD <- val[["logDet"]]
    val <- val[["factor"]]
  }
  if (corD[["M"]] > 1) {
    val <- split(val, rep(1:corD[["M"]], (corD[["len"]])^2))
    val <- lapply(val, function(el) {
      nel <- round(sqrt(length(el)))
      array(el, c(nel, nel))
    })
    names(val) <- names(corD[["len"]])
    val <- as.list(val)
  } else {
    val <- array(val, c(corD[["N"]], corD[["N"]]))
  }
  attr(val, "logDet") <- lD
  val
}

###*# Methods for standard generics

coef.corSpatial <-
  function(object, unconstrained = TRUE, ...)
{
  if (attr(object, "fixed") && unconstrained) {
    return(numeric(0))
  }
  val <- as.vector(object)
  if (length(val) == 0) {               # uninitialized
    return(val)
  }
  if (!unconstrained) {
    val <- exp(val)
    if (attr(object, "nugget")) val[2] <- val[2]/(1+val[2])
  }
  names(val) <- if(attr(object, "nugget")) c("range", "nugget") else "range"
  val
}

`coef<-.corSpatial` <- function(object, ..., value)
{
  if (length(value) != length(object)) {
    stop("cannot change the length of the parameter after initialization")
  }
  object[] <- value
  corD <- attr(object, "Dim")
  ## updating the factor list and logDet
  aux <- .C(spatial_factList,
	    as.double(as.vector(object)),
	    as.integer(attr(object, "nugget")),
	    as.double(unlist(getCovariate(object))),
	    as.integer(unlist(corD)),
	    as.double(attr(object, "minD")),
	    factor = double(corD[["sumLenSq"]]),
	    logDet = double(1))[c("factor", "logDet")]
  attr(object, "factor") <- aux[["factor"]]
  attr(object, "logDet") <- -aux[["logDet"]]
  object
}

Dim.corSpatial <- function(object, groups, ...)
{
  if (missing(groups)) return(attr(object, "Dim"))
  val <- Dim.corStruct(object, groups)
  val[["start"]] <-
    c(0, cumsum(val[["len"]] * (val[["len"]] - 1)/2)[-val[["M"]]])
  ## will use third component of Dim list for spClass
  names(val)[3] <- "spClass"
  val[[3]] <-
    match(class(object)[1], c("corSpher", "corExp", "corGaus", "corLin",
                              "corRatio"), 0)
  val
}

getCovariate.corSpatial <-
  function(object, form = formula(object), data)
{
  if (is.null(covar <- attr(object, "covariate"))) { # need to calculate it
    if (missing(data)) {
      stop("need data to calculate covariate")
    }
    covForm <- getCovariateFormula(form)
    covar <-
        if (length(all.vars(covForm)) > 0) { # covariate present
            if (attr(terms(covForm), "intercept") == 1) {
                covForm <- eval(substitute(~ CV - 1, list(CV = covForm[[2]])))
            }
            as.data.frame(unclass(
              model.matrix(covForm,
                           model.frame(covForm, data,
                                       drop.unused.levels = TRUE))))
        } ## else NULL

    covar <-
      if (!is.null(getGroupsFormula(form))) { # by groups
        grps <- getGroups(object, data = data)
        grps <-
          if (is.null(covar)) {
            lapply(split(grps, grps), function(x) as.vector(dist(seq_along(x))))
          } else {
            lapply(split(covar, grps),
                   function(el, metric) {
                     el <- as.matrix(el)
                     if (nrow(el) > 1) as.vector(dist(el, metric)) else numeric(0)
                   }, metric = attr(object, "metric"))
          }
        grps[lengths(grps) > 0]# no 1-obs groups
      } else { # no groups
        as.vector(
          if (is.null(covar))
            dist(1:nrow(data))
          else
            dist(as.matrix(covar), method = attr(object, "metric")))
      }
    if (any(unlist(covar) == 0)) {
      stop("cannot have zero distances in \"corSpatial\"")
    }
  }
  covar
}

Initialize.corSpatial <-
  function(object, data, ...)
{
  if (!is.null(attr(object, "minD"))) { #already initialized
    return(object)
  }
  object <- Initialize.corStruct(object, data)
  nug <- attr(object, "nugget")

  val <- as.vector(object)
  if (length(val) > 0) {		# initialized
    if (val[1] <= 0) {
      stop("'range' must be > 0 in \"corSpatial\" initial value")
    }
    if (nug) {				# with nugget effect
      if (length(val) == 1) {		# assuming nugget effect not given
	val <- c(val, 0.1)		# setting it to 0.1
      } else {
	if (length(val) != 2) {
	  stop("initial value for \"corSpatial\" parameters of wrong dimension")
	}
      }
      if ((val[2] <= 0) || (val[2] >= 1)) {
	stop("initial value of nugget ratio must be between 0 and 1")
      }
    } else {				# only range parameter
      if (length(val) != 1) {
	stop("initial value for \"corSpatial\" parameters of wrong dimension")
      }
    }
  } else {
    val <- min(unlist(attr(object, "covariate"))) * 0.9
    if (nug) val <- c(val, 0.1)
  }
  val[1] <- log(val[1])
  if (nug) val[2] <- log(val[2]/(1 - val[2]))
  oldAttr <- attributes(object)
  object <- val
  attributes(object) <- oldAttr
  attr(object, "minD") <- min(unlist(attr(object, "covariate")))
  attr(object, "factor") <- corFactor(object)
  attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
  object
}

recalc.corSpatial <-
  function(object, conLin, ...)
{
  val <-
    .C(spatial_recalc,
       Xy = as.double(conLin[["Xy"]]),
       as.integer(unlist(Dim(object))),
       as.integer(ncol(conLin[["Xy"]])),
       as.double(as.vector(object)),
       as.double(unlist(getCovariate(object))),
       as.double(attr(object, "minD")),
       as.integer(attr(object, "nugget")),
       logLik = double(1))[c("Xy", "logLik")]
  conLin[["Xy"]][] <- val[["Xy"]]
  conLin[["logLik"]] <- conLin[["logLik"]] + val[["logLik"]]
  conLin
}

Variogram.corSpatial <-
  function(object, distance = NULL, sig2 = 1, length.out = 50, FUN, ...)
{
  if (is.null(distance)) {
    rangeDist <- range(unlist(getCovariate(object)))
    distance <- seq(rangeDist[1], rangeDist[2], length.out = length.out)
  }
  params <- coef(object, unconstrained = FALSE)
  if (length(params) == 1) {            # no nugget effect
    rang <- params
    nugg <- 0
  } else {                              # nugget effect
    rang <- params[1]
    nugg <- params[2]
  }
  val <- data.frame(variog = sig2 * (nugg + (1 - nugg) * FUN(distance, rang)),
                    dist = distance)
  class(val) <- c("Variogram", "data.frame")
  val
}

###*# corExp - exponential spatial correlation structure

corExp <-
  ## Constructor for the corExp class
  function(value = numeric(0), form = ~ 1, nugget = FALSE,
	   metric = c("euclidean", "maximum", "manhattan"), fixed = FALSE)
{
  attr(value, "formula") <- form
  attr(value, "nugget") <- nugget
  attr(value, "metric") <- match.arg(metric)
  attr(value, "fixed") <- fixed
  class(value) <- c("corExp", "corSpatial", "corStruct")
  value
}

###*# Methods for standard generics

summary.corExp <-
  function(object, structName = "Exponential spatial correlation", ...)
{
  summary.corStruct(object, structName)
}

Variogram.corExp <-
  function(object, distance = NULL, sig2 = 1, length.out = 50, ...)
{
  Variogram.corSpatial(object, distance, sig2, length.out,
                       function(x, y) { 1 - exp(-x/y) })
}

###*# corGaus - Gaussian spatial correlation structure

corGaus <-
  ## Constructor for the corGaus class
  function(value = numeric(0), form = ~ 1, nugget = FALSE,
	   metric = c("euclidean", "maximum", "manhattan"), fixed = FALSE)
{
  attr(value, "formula") <- form
  attr(value, "nugget") <- nugget
  attr(value, "metric") <- match.arg(metric)
  attr(value, "fixed") <- fixed
  class(value) <- c("corGaus", "corSpatial", "corStruct")
  value
}

###*# Methods for standard generics

summary.corGaus <-
  function(object, structName = "Gaussian spatial correlation", ...)
{
  summary.corStruct(object, structName)
}

Variogram.corGaus <-
  function(object, distance = NULL, sig2 = 1, length.out = 50, ...)
{
  Variogram.corSpatial(object, distance, sig2, length.out,
                       function(x, y){ 1 - exp(-(x/y)^2) })
}

###*# corLin - Linear spatial correlation structure

corLin <-
  ## Constructor for the corLin class
  function(value = numeric(0), form = ~ 1, nugget = FALSE,
	   metric = c("euclidean", "maximum", "manhattan"), fixed = FALSE)
{
  attr(value, "formula") <- form
  attr(value, "nugget") <- nugget
  attr(value, "metric") <- match.arg(metric)
  attr(value, "fixed") <- fixed
  class(value) <- c("corLin", "corSpatial", "corStruct")
  value
}

###*# Methods for standard generics

coef.corLin <-
  function(object, unconstrained = TRUE, ...)
{
  val <- NextMethod()
  if (!unconstrained) val[1] <- val[1] + attr(object, "minD")
  val
}

Initialize.corLin <-
  function(object, data, ...)
{
  if (!is.null(attr(object, "minD"))) { #already initialized
    return(object)
  }
  object <- Initialize.corStruct(object, data)
  nug <- attr(object, "nugget")

  minD <- min(unlist(attr(object, "covariate")))
  val <- as.vector(object)
  if (length(val) > 0) {		# initialized
    if (val[1] <= 0) {
      stop("'range' must be > 0 in \"corLin\" initial value")
    }
    if (val[1] <= minD) {
      warning("initial value for 'range' less than minimum distance. Setting it to 1.1 * min(distance)")
      val[1] <- 1.1 * minD
    }
    if (nug) {				# with nugget effect
      if (length(val) == 1) {		# assuming nugget effect not given
	val <- c(val, 0.1)		# setting it to 0.1
      } else {
	if (length(val) != 2) {
	  stop("initial value for \"corLin\" parameters of wrong dimension")
	}
      }
      if ((val[2] <= 0) || (val[2] >= 1)) {
	stop("initial value of nugget ratio must be between 0 and 1")
      }
    } else {				# only range parameter
      if (length(val) != 1) {
	stop("initial value for \"corLin\" parameters of wrong dimension")
      }
    }
  } else {
    val <- minD * 1.1
    if (nug) val <- c(val, 0.1)
  }
  val[1] <- log(val[1] - minD)
  if (nug) val[2] <- log(val[2]/(1 - val[2]))
  oldAttr <- attributes(object)
  object <- val
  attributes(object) <- oldAttr
  attr(object, "minD") <- minD
  attr(object, "factor") <- corFactor(object)
  attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
  object
}

summary.corLin <-
  function(object, structName = "Linear spatial correlation", ...)
{
  summary.corStruct(object, structName)
}

Variogram.corLin <-
  function(object, distance = NULL, sig2 = 1, length.out = 50, ...)
{
  Variogram.corSpatial(object, distance, sig2, length.out,
                       function(x, y) { pmin(x/y, 1) })
}

###*# corRatio - rational quadratic spatial correlation structure

corRatio <-
  ## Constructor for the corRational class
  function(value = numeric(0), form = ~ 1, nugget = FALSE,
	   metric = c("euclidean", "maximum", "manhattan"), fixed = FALSE)
{
  attr(value, "formula") <- form
  attr(value, "nugget") <- nugget
  attr(value, "metric") <- match.arg(metric)
  attr(value, "fixed") <- fixed
  class(value) <- c("corRatio", "corSpatial", "corStruct")
  value
}

###*# Methods for standard generics

summary.corRatio <-
  function(object, structName = "Rational quadratic spatial correlation", ...)
{
  summary.corStruct(object, structName)
}

Variogram.corRatio <-
  function(object, distance = NULL, sig2 = 1, length.out = 50, ...)
{
  Variogram.corSpatial(object, distance, sig2, length.out,
                       function(x, y) {
                         x <- (x/y)^2
                         x/(1+x)
                       })
}

###*# corSpher - spherical spatial correlation structure

corSpher <-
  ## Constructor for the corSpher class
  function(value = numeric(0), form = ~ 1, nugget = FALSE,
	   metric = c("euclidean", "maximum", "manhattan"), fixed = FALSE)
{
  attr(value, "formula") <- form
  attr(value, "nugget") <- nugget
  attr(value, "metric") <- match.arg(metric)
  attr(value, "fixed") <- fixed
  class(value) <- c("corSpher", "corSpatial", "corStruct")
  value
}

###*# Methods for standard generics

coef.corSpher <-
  function(object, unconstrained = TRUE, ...)
{
  val <- NextMethod()
  if (!unconstrained) val[1] <- val[1] + attr(object, "minD")
  val
}

Initialize.corSpher <-
  function(object, data, ...)
{
  if (!is.null(attr(object, "minD"))) { #already initialized
    return(object)
  }
  object <- Initialize.corStruct(object, data)
  nug <- attr(object, "nugget")

  minD <- min(unlist(attr(object, "covariate")))
  val <- as.vector(object)
  if (length(val) > 0) {		# initialized
    if (val[1] <= 0) {
      stop("range must be > 0 in \"corSpher\" initial value")
    }
    if (val[1] <= minD) {
      warning("initial value for 'range' less than minimum distance. Setting it to 1.1 * min(distance)")
      val[1] <- 1.1 * minD
    }
    if (nug) {				# with nugget effect
      if (length(val) == 1) {		# assuming nugget effect not given
	val <- c(val, 0.1)		# setting it to 0.1
      } else {
	if (length(val) != 2) {
	  stop("initial value for \"corSpher\" parameters of wrong dimension")
	}
      }
      if ((val[2] <= 0) || (val[2] >= 1)) {
	stop("initial value of nugget ratio must be between 0 and 1")
      }
    } else {				# only range parameter
      if (length(val) != 1) {
	stop("initial value for \"corSpher\" parameters of wrong dimension")
      }
    }
  } else {
    val <- minD * 1.1
    if (nug) val <- c(val, 0.1)
  }
  val[1] <- log(val[1] - minD)
  if (nug) val[2] <- log(val[2]/(1 - val[2]))
  oldAttr <- attributes(object)
  object <- val
  attributes(object) <- oldAttr
  attr(object, "minD") <- minD
  attr(object, "factor") <- corFactor(object)
  attr(object, "logDet") <- -attr(attr(object, "factor"), "logDet")
  object
}

summary.corSpher <-
  function(object, structName = "Spherical spatial correlation", ...)
{
  summary.corStruct(object, structName)
}

Variogram.corSpher <-
  function(object, distance = NULL, sig2 = 1, length.out = 50, ...)
{
  Variogram.corSpatial(object, distance, sig2, length.out,
                       function(x, y) {
                         x <- pmin(x/y, 1)
                         1.5 * x - 0.5 * x^3
                       })
}

####*# corWave - Wave spatial correlation structure

#corWave <-
#  ## Constructor for the corWave class
#  function(value = numeric(0), form = ~ 1, nugget = FALSE,
#	   metric = c("euclidean", "maximum", "manhattan"))
#{
#  attr(value, "formula") <- form
#  attr(value, "nugget") <- nugget
#  attr(value, "metric") <- match.arg(metric)
#  class(value) <- c("corWave", "corSpatial", "corStruct")
#  value
#}

####*# Methods for standard generics

#summary.corWave <-
#  function(object, structName = "Wave spatial correlation")
#{
#  summary.corStruct(object, structName)
#}


##*## Beginning of epilogue
### This file is automatically placed in Outline minor mode.
### The file is structured as follows:
### Chapters:     ^L #
### Sections:    ##*##
### Subsections: ###*###
### Components:  non-comment lines flushed left
###              Random code beginning with a ####* comment

### Local variables:
### mode: outline-minor
### outline-regexp: "\^L\\|\\`#\\|##\\*\\|###\\*\\|[a-zA-Z]\\|\\\"[a-zA-Z]\\|####\\*"
### ess-indent-offset: 2
### End: