if((Rv <- getRversion()) < "3.2.1") {
    lengths <- function (x, use.names = TRUE) vapply(x, length, 1L, USE.NAMES = use.names)
    if(Rv < "3.1.0") {
        anyNA <- function(x) any(is.na(x))
        if(Rv < "3.0.0") {
            rep_len <- function(x, length.out) rep(x, length.out=length.out)
            if(Rv < "2.15")
                paste0 <- function(...) paste(..., sep = '')
        }
    }
} ## R < 3.2.1
rm(Rv)

## From Matrix package  isDiagonal(.) :
all0 <- function(x) !anyNA(x) && all(!x)
.isDiagonal.sq.matrix <- function(M, n = dim(M)[1L])
    all0(M[rep_len(c(FALSE, rep.int(TRUE,n)), n^2)])


### Utilities for parsing and manipulating mixed-model formulas

##' deparse(.) returning \bold{one} string
##' @note Protects against the possibility that results from deparse() will be
##'       split after 'width.cutoff' (by default 60, maximally 500)
safeDeparse <- function(x, collapse=" ") paste(deparse(x, 500L), collapse=collapse)
abbrDeparse <- function(x, width=60) {
    r <- deparse(x, width)
    if(length(r) > 1) paste(r[1], "...") else r
}

##' @param bars result of findbars
barnames <- function(bars) vapply(bars, function(x) safeDeparse(x[[3]]), "")

makeFac <- function(x,char.only=FALSE) {
    if (!is.factor(x) && (!char.only || is.character(x))) factor(x) else x
}

factorize <- function(x,frloc,char.only=FALSE) {
    ## convert grouping variables to factors as necessary
    ## TODO: variables that are *not* in the data frame are
    ##  not converted -- these could still break, e.g. if someone
    ##  tries to use the : operator
    ## TODO: some sensible tests for drop.unused.levels
    ##       (not actually used, but could come in handy)
    for (i in all.vars(RHSForm(x))) {
        if (!is.null(curf <- frloc[[i]]))
            frloc[[i]] <- makeFac(curf,char.only)
    }
    return(frloc)
}

##' @param x a language object of the form  effect | groupvar
##' @param frloc model frame
##' @param drop.unused.levels (logical)
##' @return list containing grouping factor, sparse model matrix, number of levels, names
mkBlist <- function(x,frloc, drop.unused.levels=TRUE) {
    frloc <- factorize(x,frloc)
    ## try to evaluate grouping factor within model frame ...
    if (is.null(ff <- tryCatch(eval(substitute(makeFac(fac),
                                               list(fac = x[[3]])), frloc),
                error=function(e) NULL)))
        stop("couldn't evaluate grouping factor ",
             deparse(x[[3]])," within model frame:",
             " try adding grouping factor to data ",
             "frame explicitly if possible",call.=FALSE)
    if (all(is.na(ff)))
        stop("Invalid grouping factor specification, ",
             deparse(x[[3]]),call.=FALSE)
    ## NB: *also* silently drops <NA> levels - and mkReTrms() and hence
    ##     predict.merMod() have relied on that property  :
    if (drop.unused.levels) ff <- factor(ff, exclude=NA)
    nl <- length(levels(ff))
    ## this section implements eq. 6 of the JSS lmer paper
    ## model matrix based on LHS of random effect term (X_i)
    ##    x[[2]] is the LHS (terms) of the a|b formula
    mm <- model.matrix(eval(substitute( ~ foo, list(foo = x[[2]]))), frloc)
    ## nc <- ncol(mm)
    ## nseq <- seq_len(nc)
    ## this is J^T (see p. 9 of JSS lmer paper)
    ## construct indicator matrix for groups by observations
    ## use fac2sparse() rather than as() to allow *not* dropping
    ## unused levels where desired
    sm <- fac2sparse(ff, to = "d",
                     drop.unused.levels = drop.unused.levels)
    ## looks like we don't have to filter NAs explicitly any more ...
    ## sm <- as(ff,"sparseMatrix")
    ## sm <- KhatriRao(sm[,!is.na(ff),drop=FALSE],t(mm[!is.na(ff),,drop=FALSE]))
    sm <- KhatriRao(sm, t(mm))
    dimnames(sm) <- list(
        rep(levels(ff),each=ncol(mm)),
        rownames(mm))
    list(ff = ff, sm = sm, nl = nl, cnms = colnames(mm))
}

##' From the result of \code{\link{findbars}} applied to a model formula and
##' and the evaluation frame, create the model matrix, etc. associated with
##' random-effects terms.  See the description of the returned value for a
##' detailed list.
##'
##' @title Create Z, Lambda, Lind, etc.
##' @param bars a list of parsed random-effects terms
##' @param fr a model frame in which to evaluate these terms
##' @return a list with components
##' \item{Zt}{transpose of the sparse model matrix for the random effects}
##' \item{Lambdat}{transpose of the sparse relative covariance factor}
##' \item{Lind}{an integer vector of indices determining the mapping of the
##'     elements of the \code{theta} to the \code{"x"} slot of \code{Lambdat}}
##' \item{theta}{initial values of the covariance parameters}
##' \item{lower}{lower bounds on the covariance parameters}
##' \item{flist}{list of grouping factors used in the random-effects terms}
##' \item{cnms}{a list of column names of the random effects according to
##'     the grouping factors}
##' @importFrom Matrix sparseMatrix drop0
##' @importMethodsFrom Matrix coerce rbind
##' @family utilities
##' @export
mkReTrms <- function(bars, fr, drop.unused.levels=TRUE) {
  if (!length(bars))
    stop("No random effects terms specified in formula",call.=FALSE)
  stopifnot(is.list(bars), vapply(bars, is.language, NA),
            inherits(fr, "data.frame"))
  names(bars) <- barnames(bars)
  term.names <- vapply(bars, safeDeparse, "")
  ## get component blocks
  blist <- lapply(bars, mkBlist, fr, drop.unused.levels)
  nl <- vapply(blist, `[[`, 0L, "nl")   # no. of levels per term
                                        # (in lmer jss:  \ell_i)

  ## order terms stably by decreasing number of levels in the factor
  if (any(diff(nl) > 0)) {
    ord <- rev(order(nl))
    blist      <- blist     [ord]
    nl         <- nl        [ord]
    term.names <- term.names[ord]
  }
  Ztlist <- lapply(blist, `[[`, "sm")
  Zt <- do.call(rbind, Ztlist)  ## eq. 7, JSS lmer paper
  names(Ztlist) <- term.names
  q <- nrow(Zt)

  ## Create and install Lambdat, Lind, etc.  This must be done after
  ## any potential reordering of the terms.
  cnms <- lapply(blist, `[[`, "cnms")   # list of column names of the
                                        # model matrix per term
  nc <- lengths(cnms)	                # no. of columns per term
                                        # (in lmer jss:  p_i)
  nth <- as.integer((nc * (nc+1))/2)	# no. of parameters per term
                                        # (in lmer jss:  ??)
  nb <- nc * nl			        # no. of random effects per term
                                        # (in lmer jss:  q_i)
  ## eq. 5, JSS lmer paper
  if (sum(nb) != q) {
      stop(sprintf("total number of RE (%d) not equal to nrow(Zt) (%d)",
                   sum(nb),q))
  }
  boff <- cumsum(c(0L, nb))		# offsets into b
  thoff <- cumsum(c(0L, nth))		# offsets into theta
  ### FIXME: should this be done with cBind and avoid the transpose
  ### operator?  In other words should Lambdat be generated directly
  ### instead of generating Lambda first then transposing?
  Lambdat <-
    t(do.call(sparseMatrix,
              do.call(rbind,
                      lapply(seq_along(blist), function(i)
                      {
                        mm <- matrix(seq_len(nb[i]), ncol = nc[i],
                                     byrow = TRUE)
                        dd <- diag(nc[i])
                        ltri <- lower.tri(dd, diag = TRUE)
                        ii <- row(dd)[ltri]
                        jj <- col(dd)[ltri]
                        ## unused: dd[cbind(ii, jj)] <- seq_along(ii)
                        data.frame(i = as.vector(mm[, ii]) + boff[i],
                                   j = as.vector(mm[, jj]) + boff[i],
                                   x = as.double(rep.int(seq_along(ii),
                                                         rep.int(nl[i], length(ii))) +
                                                   thoff[i]))
                      }))))
  thet <- numeric(sum(nth))
  ll <- list(Zt = drop0(Zt), theta = thet, Lind = as.integer(Lambdat@x),
             Gp = unname(c(0L, cumsum(nb))))
  ## lower bounds on theta elements are 0 if on diagonal, else -Inf
  ll$lower <- -Inf * (thet + 1)
  ll$lower[unique(diag(Lambdat))] <- 0
  ll$theta[] <- is.finite(ll$lower) # initial values of theta are 0 off-diagonal, 1 on
  Lambdat@x[] <- ll$theta[ll$Lind]  # initialize elements of Lambdat
  ll$Lambdat <- Lambdat
  # massage the factor list
  fl <- lapply(blist, `[[`, "ff")
  # check for repeated factors
  fnms <- names(fl)
  if (length(fnms) > length(ufn <- unique(fnms))) {
    fl <- fl[match(ufn, fnms)]
    asgn <- match(fnms, ufn)
  } else asgn <- seq_along(fl)
  names(fl) <- ufn
  ## DON'T need fl to be a data.frame ...
  ## fl <- do.call(data.frame, c(fl, check.names = FALSE))
  attr(fl, "assign") <- asgn
  ll$flist <- fl
  ll$cnms <- cnms
  ll$Ztlist <- Ztlist
  ll
} ## {mkReTrms}

##' Create an lmerResp, glmResp or nlsResp instance
##'
##' @title Create an lmerResp, glmResp or nlsResp instance
##' @param fr a model frame
##' @param REML logical scalar, value of REML for an lmerResp instance
##' @param family the optional glm family (glmResp only)
##' @param nlenv the nonlinear model evaluation environment (nlsResp only)
##' @param nlmod the nonlinear model function (nlsResp only)
##' @param ... where to look for response information if \code{fr} is missing.
##'   Can contain a model response, \code{y}, offset, \code{offset}, and weights,
##'   \code{weights}.
##' @return an lmerResp or glmResp or nlsResp instance
##' @family utilities
##' @export
mkRespMod <- function(fr, REML=NULL, family = NULL, nlenv = NULL, nlmod = NULL, ...)
{
    if(!missing(fr)) {
        y <- model.response(fr)
        offset <- model.offset(fr)
        weights <- model.weights(fr)
        N <- n <- nrow(fr)
        etastart_update <- model.extract(fr, "etastart")
        mustart_update <- model.extract(fr, "mustart")
    } else {
        fr <- list(...)
        y <- fr$y
        N <- n <- NROW(y)
        offset <- fr$offset
        weights <- fr$weights
        etastart_update <- fr$etastart
        mustart_update <- fr$mustart
    }
    if(length(dim(y)) == 1L)
	y <- drop(y) ## avoid problems with 1D arrays and keep names

    if(isGLMM <- !is.null(family))
        stopifnot(inherits(family, "family"))
   ## FIXME: may need to add X, or pass it somehow, if we want to use glm.fit

    ## test for non-numeric response here to avoid later
    ## confusing error messages from deeper machinery
    if (!is.null(y)) { ## 'y' may be NULL if we're doing simulation
	if(!(is.numeric(y) ||
	    ((is.binom <- isGLMM && family$family == "binomial") &&
                (is.factor(y) || is.logical(y))))) {
	    if (is.binom)
		stop("response must be numeric or factor")
	    else {
                if (is.logical(y))
                    y <- as.integer(y)
                else stop("response must be numeric")
            }
	}
	if(!all(is.finite(y)))
	    stop("NA/NaN/Inf in 'y'") # same msg as from lm.fit()
    }

    rho <- new.env()
    rho$y <- if (is.null(y)) numeric(0) else y
    if (!is.null(REML)) rho$REML <- REML
    rho$etastart <- etastart_update
    rho$mustart <- mustart_update
    rho$start <- attr(fr,"start")
    if (!is.null(nlenv)) {
        stopifnot(is.language(nlmod),
                  is.environment(nlenv),
                  is.numeric(val <- eval(nlmod, nlenv)),
                  length(val) == n,
		  ## FIXME?  Restriction, not present in ole' nlme():
                  is.matrix(gr <- attr(val, "gradient")),
                  is.numeric(gr),
                  nrow(gr) == n,
                  !is.null(pnames <- colnames(gr)))
        N <- length(gr)
        rho$mu <- as.vector(val)
        rho$sqrtXwt <- as.vector(gr)
        rho$gam <- ## FIXME more efficient  mget(pnames, envir=nlenv)
            unname(unlist(lapply(pnames,
                                 function(nm) get(nm, envir=nlenv))))
    }
    rho$offset <- if (!is.null(offset)) {
        if (length(offset) == 1L) offset <- rep.int(offset, N)
        else stopifnot(length(offset) == N)
        unname(offset)
    } else rep.int(0, N)
    rho$weights <- if (!is.null(weights)) {
        stopifnot(length(weights) == n, all(weights >= 0))
        unname(weights)
    } else rep.int(1, n)

    if(isGLMM) {
        ## need weights for initializing evaluation
        rho$nobs <- n
        ## allow trivial objects, e.g. for simulation
        if (length(y)>0) eval(family$initialize, rho)
        ## ugh. this *is* necessary;
        ##  family$initialize *ignores* mustart in env, overwrites!
        ## see ll 180-182 of src/library/stats/R/glm.R
        ## https://github.com/wch/r-source/search?utf8=%E2%9C%93&q=mukeep
        if (!is.null(mustart_update)) rho$mustart <- mustart_update
        ## family$initialize <- NULL     # remove clutter from str output
        ll <- as.list(rho)
        ans <- do.call(new, c(list(Class="glmResp", family=family),
                              ll[setdiff(names(ll), c("m", "nobs", "mustart"))]))
        if (length(y)>0)
            ans$updateMu(if (!is.null(es <- etastart_update)) es else                                       family$linkfun(rho$mustart))
        ans
    } else if (is.null(nlenv)) ## lmer
        do.call(lmerResp$new, as.list(rho))
    else ## nlmer
        do.call(nlsResp$new,
                c(list(nlenv=nlenv,
                       nlmod=substitute(~foo, list(foo=nlmod)),
                       pnames=pnames), as.list(rho)))
}

##' From the right hand side of a formula for a mixed-effects model,
##' determine the pairs of expressions that are separated by the
##' vertical bar operator.  Also expand the slash operator in grouping
##' factor expressions and expand terms with the double vertical bar operator
##' into separate, independent random effect terms.
##'
##' @title Determine random-effects expressions from a formula
##' @seealso \code{\link{formula}}, \code{\link{model.frame}}, \code{\link{model.matrix}}.
##' @param term a mixed-model formula
##' @return pairs of expressions that were separated by vertical bars
##' @section Note: This function is called recursively on individual
##' terms in the model, which is why the argument is called \code{term} and not
##' a name like \code{form}, indicating a formula.
##' @example
##' findbars(f1 <- Reaction ~ Days + (Days|Subject))
##' ## => list( Days | Subject )
##' findbars(y ~ Days + (1|Subject) + (0+Days|Subject))
##' ## => list of length 2:  list ( 1 | Subject ,  0+Days|Subject)
##' findbars(~ 1 + (1|batch/cask))
##' ## => list of length 2:  list ( 1 | cask:batch ,  1 | batch)
##' identical(findbars(~ 1 + (Days || Subject)),
##'     findbars(~ 1 + (1|Subject) + (0+Days|Subject)))
##' \dontshow{
##' stopifnot(identical(findbars(f1),
##'                     list(expression(Days | Subject)[[1]])))
##' }
##' @family utilities
##' @keywords models utilities
##' @export
findbars <- function(term)
{
    ## Recursive function applied to individual terms
    fb <- function(term)
    {
	if (is.name(term) || !is.language(term)) return(NULL)
	if (term[[1]] == as.name("(")) return(fb(term[[2]]))
	stopifnot(is.call(term))
	if (term[[1]] == as.name('|')) return(term)
	if (length(term) == 2) return(fb(term[[2]]))
	c(fb(term[[2]]), fb(term[[3]]))
    }
    ## Expand any slashes in the grouping factors returned by fb
    expandSlash <- function(bb)
    {
	## Create the interaction terms for nested effects
	makeInteraction <- function(x)
	{
	    if (length(x) < 2) return(x)
	    trm1 <- makeInteraction(x[[1]])
	    trm11 <- if(is.list(trm1)) trm1[[1]] else trm1
	    list(substitute(foo:bar, list(foo=x[[2]], bar = trm11)), trm1)
	}
	## Return the list of '/'-separated terms
	slashTerms <- function(x)
	{
	    if (!("/" %in% all.names(x))) return(x)
	    if (x[[1]] != as.name("/"))
		stop("unparseable formula for grouping factor",call.=FALSE)
	    list(slashTerms(x[[2]]), slashTerms(x[[3]]))
	}

	if (!is.list(bb))
	    expandSlash(list(bb))
	else
	    unlist(lapply(bb, function(x) {
		if (length(x) > 2 && is.list(trms <- slashTerms(x[[3]])))
		    ## lapply(unlist(...)) - unlist returns a flattened list
		    lapply(unlist(makeInteraction(trms)),
			   function(trm) substitute(foo|bar, list(foo = x[[2]], bar = trm)))
		else x
	    }))
    }## {expandSlash}

    modterm <- expandDoubleVerts(
	if(is(term, "formula")) term[[length(term)]] else term)
    expandSlash(fb(modterm))
}

##' From the right hand side of a formula for a mixed-effects model,
##' expand terms with the double vertical bar operator
##' into separate, independent random effect terms.
##'
##' @title Expand terms with \code{'||'} notation into separate \code{'|'} terms
##' @seealso \code{\link{formula}}, \code{\link{model.frame}}, \code{\link{model.matrix}}.
##' @param term a mixed-model formula
##' @return the modified term
##' @family utilities
##' @keywords models utilities
##' @export
expandDoubleVerts <- function(term)
{
    expandDoubleVert <- function(term) {
        frml <- formula(substitute(~x,list(x=term[[2]])))
        ## FIXME: do this without paste and deparse if possible!
	## need term.labels not all.vars to capture interactions too:
	newtrms <- paste0("0+", attr(terms(frml), "term.labels"))
	if(attr(terms(frml), "intercept")!=0)
	    newtrms <- c("1", newtrms)

	as.formula(paste("~(",
			 paste(vapply(newtrms, function(trm)
				      paste0(trm, "|", deparse(term[[3]])), ""),
			       collapse=")+("), ")"))[[2]]
    }

    if (!is.name(term) && is.language(term)) {
	if (term[[1]] == as.name("(")) {
	    term[[2]] <- expandDoubleVerts(term[[2]])
	}
	stopifnot(is.call(term))
	if (term[[1]] == as.name('||'))
	    return( expandDoubleVert(term) )
	## else :
	term[[2]] <- expandDoubleVerts(term[[2]])
	if (length(term) != 2) {
	    if(length(term) == 3)
		term[[3]] <- expandDoubleVerts(term[[3]])
	}
    }
    term
}



##' Remove the random-effects terms from a mixed-effects formula,
##' thereby producing the fixed-effects formula.
##'
##' @title Omit terms separated by vertical bars in a formula
##' @param term the right-hand side of a mixed-model formula
##' @return the fixed-effects part of the formula
##' @section Note: This function is called recursively on individual
##' terms in the model, which is why the argument is called \code{term} and not
##' a name like \code{form}, indicating a formula.
##' @examples
##' nobars(Reaction ~ Days + (Days|Subject)) ## => Reaction ~ Days
##' @seealso \code{\link{formula}}, \code{\link{model.frame}}, \code{\link{model.matrix}}.
##' @family utilities
##' @keywords models utilities
##' @export
nobars <- function(term) {
    nb <- nobars_(term)  ## call recursive version
    if (is(term,"formula") && length(term)==3 && is.symbol(nb)) {
        ## called with two-sided RE-only formula:
        ##    construct response~1 formula
        nb <- reformulate("1",response=deparse(nb))
    }
    ## called with one-sided RE-only formula, or RHS alone
    if (is.null(nb)) {
        nb <- if (is(term,"formula")) ~1 else 1
    }
    nb
}

nobars_ <- function(term)
{
    if (!anyBars(term)) return(term)
    if (isBar(term)) return(NULL)
    if (isAnyArgBar(term)) return(NULL)
    if (length(term) == 2) {
        nb <- nobars_(term[[2]])
        if(is.null(nb)) return(NULL)
        term[[2]] <- nb
        return(term)
    }
    nb2 <- nobars_(term[[2]])
    nb3 <- nobars_(term[[3]])
    if (is.null(nb2)) return(nb3)
    if (is.null(nb3)) return(nb2)
    term[[2]] <- nb2
    term[[3]] <- nb3
    term
}

isBar <- function(term) {
    if(is.call(term)) {
        if((term[[1]] == as.name("|")) || (term[[1]] == as.name("||"))) {
            return(TRUE)
        }
    }
    FALSE
}

isAnyArgBar <- function(term) {
    if ((term[[1]] != as.name("~")) && (term[[1]] != as.name("("))) {
        for(i in seq_along(term)) {
            if(isBar(term[[i]])) return(TRUE)
        }
    }
    FALSE
}

anyBars <- function(term) {
    any(c('|','||') %in% all.names(term))
}


##' Substitute the '+' function for the '|' and '||' function in a mixed-model
##' formula.  This provides a formula suitable for the current
##' model.frame function.
##'
##' @title "Sub[stitute] Bars"
##' @param term a mixed-model formula
##' @return the formula with all |  and || operators replaced by +
##' @section Note: This function is called recursively on individual
##' terms in the model, which is why the argument is called \code{term} and not
##' a name like \code{form}, indicating a formula.
##' @examples
##' subbars(Reaction ~ Days + (Days|Subject)) ## => Reaction ~ Days + (Days + Subject)
##' @seealso \code{\link{formula}}, \code{\link{model.frame}}, \code{\link{model.matrix}}.
##' @family utilities
##' @keywords models utilities
##' @export
subbars <- function(term)
{
    if (is.name(term) || !is.language(term)) return(term)
    if (length(term) == 2) {
        term[[2]] <- subbars(term[[2]])
        return(term)
    }
    stopifnot(length(term) >= 3)
    if (is.call(term) && term[[1]] == as.name('|'))
        term[[1]] <- as.name('+')
    if (is.call(term) && term[[1]] == as.name('||'))
        term[[1]] <- as.name('+')
    for (j in 2:length(term)) term[[j]] <- subbars(term[[j]])
    term
}

##' Does every level of f1 occur in conjunction with exactly one level
##' of f2? The function is based on converting a triplet sparse matrix
##' to a compressed column-oriented form in which the nesting can be
##' quickly evaluated.
##'
##' @title Is f1 nested within f2?
##'
##' @param f1 factor 1
##' @param f2 factor 2
##'
##' @return TRUE if factor 1 is nested within factor 2
##' @examples
##' with(Pastes, isNested(cask, batch))   ## => FALSE
##' with(Pastes, isNested(sample, batch))  ## => TRUE
##' @export
isNested <- function(f1, f2)
{
    f1 <- as.factor(f1)
    f2 <- as.factor(f2)
    stopifnot(length(f1) == length(f2))
    k <- length(levels(f1))
    sm <- as(new("ngTMatrix",
                 i = as.integer(f2) - 1L,
                 j = as.integer(f1) - 1L,
                 Dim = c(length(levels(f2)), k)),
             "CsparseMatrix")
    all(sm@p[2:(k+1L)] - sm@p[1:k] <= 1L)
}

subnms <- function(form, nms) {
    ## Recursive function applied to individual terms
    sbnm <- function(term)
    {
        if (is.name(term)) {
	    if (any(term == nms)) 0 else term
	} else switch(length(term),
	       term, ## 1
	   {   ## 2
	       term[[2]] <- sbnm(term[[2]])
	       term
	   },
	   {   ## 3
	       term[[2]] <- sbnm(term[[2]])
	       term[[3]] <- sbnm(term[[3]])
	       term
	   })
    }
    sbnm(form)
}

##' Check for a constant term (a literal 1) in an expression
##
##' In the mixed-effects part of a nonlinear model formula, a constant
##' term is not meaningful because every term must be relative to a
##' nonlinear model parameter.  This function recursively checks the
##' expressions in the formula for a a constant, calling stop() if
##' such a term is encountered.
##' @title Check for constant terms.
##' @param expr an expression
##' @return NULL.  The function is executed for its side effect.
chck1 <- function(expr) {
    if ((le <- length(expr)) == 1) {
        if (is.numeric(expr) && expr == 1)
            stop("1 is not meaningful in a nonlinear model formula")
        return()
    } else
        for (j in seq_len(le)[-1]) Recall(expr[[j]])
}

## ---> ../man/nlformula.Rd --- Manipulate a nonlinear model formula
##' @param mc matched call from the caller, with arguments 'formula','start',...
##' @return a list with components "respMod", "frame", "X", "reTrms"
nlformula <- function(mc) {
  start <- eval(mc$start, parent.frame(2L))
  if (is.numeric(start)) start <- list(nlpars = start)
  stopifnot(is.numeric(nlpars <- start$nlpars),
            lengths(nlpars) == 1L,
            length(pnames <- names(nlpars)) == length(nlpars),
            length(form <- as.formula(mc$formula)) == 3L,
            is(nlform <- eval(form[[2]]), "formula"),
            pnames %in% all.vars(nlmod <-
                as.call(nlform[[lnl <- length(nlform)]])))

  ## MM{FIXME}: fortune(106) even twice in here!
    nlform[[lnl]] <- parse(text= paste(setdiff(all.vars(form), pnames), collapse=' + '))[[1]]
    nlform <- eval(nlform)
    environment(nlform) <- environment(form)
    m <- match(c("data", "subset", "weights", "na.action", "offset"),
               names(mc), 0)
    mc <- mc[c(1, m)]
    mc$drop.unused.levels <- TRUE
    mc[[1L]] <- quote(stats::model.frame)
    mc$formula <- nlform
    fr <- eval(mc, parent.frame(2L))
    n <- nrow(fr)
    nlenv <- list2env(fr, parent=parent.frame(2L))
    lapply(pnames, function(nm) nlenv[[nm]] <- rep.int(nlpars[[nm]], n))
    respMod <- mkRespMod(fr, nlenv=nlenv, nlmod=nlmod)

    chck1(meform <- form[[3L]])
    pnameexpr <- parse(text=paste(pnames, collapse='+'))[[1]]
    nb <- nobars_(meform)  ## call ORIGINAL recursive form
    fe <- eval(substitute(~ 0 + nb + pnameexpr))
    environment(fe) <- environment(form)
    frE <- do.call(rbind, lapply(seq_along(nlpars), function(i) fr)) # rbind s copies of the frame
    for (nm in pnames) # convert these variables in fr to indicators
        frE[[nm]] <- as.numeric(rep(nm == pnames, each = n))
    X <- model.matrix(fe, frE)
    rownames(X) <- NULL

    reTrms <- mkReTrms(lapply(findbars(meform),
                              function(expr) {
                                  expr[[2]] <- substitute(0+foo, list(foo=expr[[2]]))
                                  expr
                              }), frE)
    list(respMod=respMod, frame=fr, X=X, reTrms=reTrms, pnames=pnames)
} ## {nlformula}

################################################################################
## Beginning to think about exposing tools to create devcomp lists.
## Could be useful when extending merMod objects.  Commenting them out
## however, because R CMD check is complaining:
## https://github.com/lme4/lme4/commit/8d71e439758999ea8f90eb4752487e189407ef33#commitcomment-8773017
################################################################################
##
## .dims <- function(pp, resp, nAGQ,
##                   reTrms, n, p, rcl,
##                   compDev = NULL) {
##     if(missing(rcl)) rcl <- class(resp)
##     if(missing(n)) n <- nrow(pp$V)
##     if(missing(p)) p <- ncol(pp$V)
##     c(N=nrow(pp$X), n=n, p=p, nmp=n-p,
##       nth=length(pp$theta), q=nrow(pp$Zt),
##       nAGQ=rho$nAGQ,
##       compDev=rho$compDev,
##       ## 'use scale' in the sense of whether dispersion parameter should
##       ##  be reported/used (*not* whether theta should be scaled by sigma)
##       useSc=(rcl != "glmResp" ||
##              !resp$family$family %in% c("poisson","binomial")),
##       reTrms=length(reTrms$cnms),
##       spFe=0L,
##       REML=if (rcl=="lmerResp") resp$REML else 0L,
##       GLMM=(rcl=="glmResp"),
##       NLMM=(rcl=="nlsResp"))
## }
##
## .cmp <- function(pp, resp, dims, fval,
##                  wrss, sqrLenU, pwrss,
##                  sigmaML, rcl, fac,
##                  tolPwrss = NULL,
##                  trivial.y = FALSE) {
##     if(missing(rcl)) rcl <- class(resp)
##     if(missing(fac)) fac <- as.numeric(rcl != "nlsResp")
##     if(missing(wrss)) wrss <- resp$wrss()
##     if(missing(sqrLenU)) sqrLenU <- pp$sqrL(fac)
##     if(missing(pwrss)) pwrss <- wrss + sqrLenU
##     if(missing(sigmaML)) sigmaML <- pwrss/dims['n']
##     c(ldL2=pp$ldL2(), ldRX2=pp$ldRX2(), wrss=wrss,
##       ussq=sqrLenU, pwrss=pwrss,
##       drsum=if (rcl=="glmResp" && !trivial.y) resp$resDev() else NA,
##       REML=if (rcl=="lmerResp" && resp$REML != 0L && !trivial.y)
##       opt$fval else NA,
##       ## FIXME: construct 'REML deviance' here?
##       dev=if (rcl=="lmerResp" && resp$REML != 0L || trivial.y) NA else opt$fval,
##       sigmaML=sqrt(unname(if (!dims["useSc"] || trivial.y) NA else sigmaML)),
##       sigmaREML=sqrt(unname(if (rcl!="lmerResp" || trivial.y) NA else sigmaML*(dims['n']/dims['nmp']))),
##       tolPwrss=rho$tolPwrss)
## }
################################################################################

.minimalOptinfo <- function()
    list(conv = list(opt = 0L,
                     lme4 = list(messages = character(0))))

getConv <- function(x) {
    if (!is.null(x[["conv"]])) {
        x[["conv"]]
    } else x[["convergence"]]
}

.optinfo <- function(opt, lme4conv=NULL)
    list(optimizer = attr(opt, "optimizer"),
	 control   = attr(opt, "control"),
	 derivs    = attr(opt, "derivs"),
	 conv      = list(opt = getConv(opt), lme4 = lme4conv),
	 feval     = if (is.null(opt$feval)) NA else opt$feval,
	 warnings  = attr(opt, "warnings"),
	 val       = opt$par)

##' Potentically needed in more than one place, be sure to keep consistency!
##' hack (NB families have weird names) from @aosmith16; then corrected
isNBfamily <- function(familyString)
    grepl("^Negative ?Binomial", familyString, ignore.case=TRUE)
normalizeFamilyName <- function(family) { # such as  object@resp$family
    if(isNBfamily(family$family))
	family$family <- "negative.binomial"
    family
}

##' Is it a family with no scale parameter
hasNoScale <- function(family)
    any(substr(family$family, 1L, 12L)
        == c("poisson", "binomial", "negative.bin", "Negative Bin"))



##--> ../man/mkMerMod.Rd ---Create a merMod object
##' @param rho the environment of the objective function
##' @param opt the value returned by the optimizer
##' @param reTrms reTrms list from the calling function
mkMerMod <- function(rho, opt, reTrms, fr, mc, lme4conv=NULL) {
    if(missing(mc)) mc <- match.call()
    stopifnot(is.environment(rho),
              is(pp <- rho$pp, "merPredD"),
              is(resp <- rho$resp, "lmResp"),
              is.list(opt), "par" %in% names(opt),
              c("conv", "fval") %in% substr(names(opt),1,4), ## "conv[ergence]", "fval[ues]"
              is.list(reTrms), c("flist", "cnms", "Gp", "lower") %in% names(reTrms),
              length(rcl <- class(resp)) == 1)
    n    <- nrow(pp$V)
    p    <- ncol(pp$V)
    isGLMM <- (rcl == "glmResp")
    dims <- c(N = nrow(pp$X), n=n, p=p, nmp = n-p, q = nrow(pp$Zt),
              nth = length(pp$theta),
              nAGQ= rho$nAGQ,
              compDev=rho$compDev,
              ## 'use scale' in the sense of whether dispersion parameter should
              ##  be reported/used (*not* whether theta should be scaled by sigma)
              useSc = !(isGLMM && hasNoScale(resp$family)),
              reTrms=length(reTrms$cnms),
              spFe= 0L,
              REML = if (rcl=="lmerResp") resp$REML else 0L,
              GLMM= isGLMM,
              NLMM= (rcl=="nlsResp"))
    storage.mode(dims) <- "integer"
    fac     <- as.numeric(rcl != "nlsResp")
    if (trivial.y <- (length(resp$y)==0)) {
        ## trivial model
        sqrLenU <- wrss <- pwrss <- NA
    } else {
        sqrLenU <- pp$sqrL(fac)
        wrss    <- resp$wrss()
        pwrss   <- wrss + sqrLenU
    }
    ## weights <- resp$weights
    beta    <- pp$beta(fac)
    ## rescale
    if (!is.null(sc <- attr(pp$X, "scaled:scale"))) {
        warning("auto(un)scaling not yet finished/tested")
        ## FIXME: test/handle no-intercept models
        ##   (only need to worry if we do centering as well as scaling)
        ## FIXME: adjust Hessian/vcov
        ## FIXME: where else will these changes propagate?
        ##        profiling?
        beta2 <- beta
        beta2[names(sc)] <- sc*beta2[names(sc)]
        beta <- beta2
    }
    if (!is.null(attr(pp$X, "scaled:center"))) {
        warning("auto(un)centering not yet implemented")
    }
    #sigmaML <- pwrss/sum(weights)
    sigmaML <- pwrss/n
    if (rcl != "lmerResp") {
        pars <- opt$par
        if (length(pars) > length(pp$theta)) beta <- pars[-(seq_along(pp$theta))]
    }
    cmp <- c(ldL2=pp$ldL2(), ldRX2=pp$ldRX2(), wrss=wrss,
             ussq=sqrLenU, pwrss=pwrss,
             drsum=if (rcl=="glmResp" && !trivial.y) resp$resDev() else NA,
             REML=if (rcl=="lmerResp" && resp$REML != 0L && !trivial.y)
                  opt$fval else NA,
             ## FIXME: construct 'REML deviance' here?
             dev=if (rcl=="lmerResp" && resp$REML != 0L || trivial.y) NA else opt$fval,
             sigmaML=sqrt(unname(if (!dims["useSc"] || trivial.y) NA else sigmaML)),
             sigmaREML=sqrt(unname(if (rcl!="lmerResp" || trivial.y) NA else
                                   sigmaML*(dims['n']/dims['nmp']))),
             tolPwrss=rho$tolPwrss)
    ## TODO:  improve this hack to get something in frame slot (maybe need weights, etc...)
    if(missing(fr)) fr <- data.frame(resp$y)
    new(switch(rcl, lmerResp = "lmerMod", glmResp = "glmerMod", nlsResp = "nlmerMod"),
        call=mc, frame=fr, flist=reTrms$flist, cnms=reTrms$cnms,
        Gp=reTrms$Gp, theta=pp$theta, beta=beta,
        u=if (trivial.y) rep(NA_real_,nrow(pp$Zt)) else pp$u(fac),
        lower=reTrms$lower, devcomp=list(cmp=cmp, dims=dims),
        pp=pp, resp=resp,
	optinfo = .optinfo(opt, lme4conv))
}## {mkMerMod}

## generic argument checking
## 'type': name of calling function ("glmer", "lmer", "nlmer")
##
checkArgs <- function(type,...) {
    l... <- list(...)
    if (isTRUE(l...[["sparseX"]])) warning("sparseX = TRUE has no effect at present",call.=FALSE)
    ## '...' handling up front, safe-guarding against typos ("familiy") :
    if(length(l... <- list(...))) {
        if (!is.null(l...[["family"]])) {  # call glmer if family specified
            ## we will only get here if 'family' is *not* in the arg list
            warning("calling lmer with family() is deprecated: please use glmer() instead",call.=FALSE)
            type <- "glmer"
        }
        ## Check for method argument which is no longer used
        ## (different meanings/hints depending on glmer vs lmer)
	if (!is.null(l...[["method"]])) {
            msg <- paste("Argument", sQuote("method"), "is deprecated.")
	    if (type == "lmer")
		msg <- paste(msg, "Use the REML argument to specify ML or REML estimation.")
	    else if (type == "glmer")
		msg <- paste(msg, "Use the nAGQ argument to specify Laplace (nAGQ=1) or adaptive",
			     "Gauss-Hermite quadrature (nAGQ>1).  PQL is no longer available.")
            warning(msg,call.=FALSE)
            l... <- l...[names(l...) != "method"]
        }
        if(length(l...)) {
            warning("extra argument(s) ",
                    paste(sQuote(names(l...)), collapse=", "),
                    " disregarded",call.=FALSE)
        }
    }
}

## check formula and data: return an environment suitable for evaluating
##  the formula.
## (1) if data is specified, return it
## (2) otherwise, if formula has an environment, use it
## (3) otherwise [e.g. if formula was passed as a string], try to use parent.frame(2)

## if #3 is true *and* the user is doing something tricky with nested functions,
## this may fail ...

## from http://stackoverflow.com/questions/14945274/determine-whether-evaluation-of-an-argument-will-fail-due-to-non-existence
missDataFun <- function(d) {
    ## test here
    ff <- sys.frames()
    ex <- substitute(d)
    ii <- rev(seq_along(ff))
    foundAnon <- FALSE
    for(i in ii) {
        ## trying to find where the original symbol is defined ...
        ex <- eval(substitute(substitute(x, env=sys.frames()[[n]]),
                              env = list(x = ex, n=i)))
        if (is.symbol(ex) && grepl("^\\.\\.[0-9]+$",safeDeparse(ex))) {
            ## testing for the dreaded "..1", "..2" pattern; this means
            ## we are stuck in an anonymous function somewhere ...
            ## won't be able to see whether data exist or not,
            ## but we should not flag them as missing -- definitely
            ## causes false positives
            ## might be better to check against quote(..1), quote(..2),
            ## ... but ugh ...
            foundAnon <- TRUE
            break
        }
    }
    return(!foundAnon && is.symbol(ex) && !exists(deparse(ex)))
}

checkFormulaData <- function(formula, data, checkLHS=TRUE,
                             checkData=TRUE, debug=FALSE) {
    nonexist.data <- missDataFun(data)
    wd <- tryCatch(eval(data), error = identity)
    if (wrong.data <- inherits(wd,"simpleError")) {
        wrong.data.msg <- wd$message
    }
    bad.data <- nonexist.data || wrong.data

    ## data not found (this *should* only happen with garbage input,
    ## OR when strings used as formulae -> drop1/update/etc.)
    ##
    if (bad.data || debug) {
        varex <- function(v, env) exists(v, envir=env, inherits=FALSE)
        allvars <- all.vars(as.formula(formula))
        allvarex <- function(env, vvec=allvars) all(vapply(vvec, varex, NA, env))
    }
    if (bad.data) { ## Choose helpful error message:
        if (allvarex(environment(formula))) {
            stop("bad 'data', but variables found in environment of formula: ",
                 "try specifying 'formula' as a formula rather ",
                 "than a string in the original model",call.=FALSE)
        } else {
            if (nonexist.data) {
                stop("'data' not found, and some variables missing from formula environment",call.=FALSE)
            } else {
                stop("bad 'data': ",wrong.data.msg)
            }
        }
    } else {
        denv <- ## The data as environment
            if (is.null(data)) {
                if (!is.null(ee <- environment(formula))) {
                    ee ## use environment of formula
            } else {
                ## e.g. no environment, e.g. because formula is a character vector
                ## parent.frame(2L) works because [g]lFormula (our calling environment)
                ## has been called within [g]lmer with env=parent.frame(1L)
                ## If you call checkFormulaData in some other bizarre way such that
                ## parent.frame(2L) is *not* OK, you deserve what you get
                ## calling checkFormulaData directly from the global
                ## environment should be OK, since trying to go up beyond the global
                ## environment keeps bringing you back to the global environment ...
                parent.frame(2L)
            }
        } else ## data specified
            list2env(data)
    }
    ## FIXME: set enclosing environment of denv to environment(formula), or parent.frame(2L) ?
    if (debug) {
        cat("Debugging parent frames in checkFormulaData:\n")
        ## find global environment -- could do this with sys.nframe() ?
        glEnv <- 1L
        while (!identical(parent.frame(glEnv),.GlobalEnv)) {
            glEnv <- glEnv+1L
        }
        ## where are vars?
        for (i in 1:glEnv) {
            OK <- allvarex(parent.frame(i))
            cat("vars exist in parent frame ", i)
            if (i == glEnv) cat(" (global)")
            cat(" ",OK, "\n")
        }
        cat("vars exist in env of formula ", allvarex(denv), "\n")
    } ## if (debug)

    stopifnot(!checkLHS || length(as.formula(formula,env=denv)) == 3)  ## check for two-sided formula
    return(denv)
}

## checkFormulaData <- function(formula,data) {
##     ee <- environment(formula)
##     if (is.null(ee)) {
##         ee <- parent.frame(2)
##     }
##     if (missing(data)) data <- ee
##     stopifnot(length(as.formula(formula,env=as.environment(data))) == 3)
##     return(data)
## }


##' Not exported; for tests (and examples) that can be slow;
##' Use   if(lme4:::testLevel() >= 1.) .....  see ../README.md
testLevel <- function()
    if(nzchar(s <- Sys.getenv("LME4_TEST_LEVEL")) &&
       is.finite(s <- as.numeric(s))) s else 1

##' General conditional variance-covariance matrix
##'
##' Experimental function for estimating the variance-covariance
##' matrix of the random effects, conditional on the observed data
##' and at the (RE)ML estimate of the fixed effects and covariance
##' parameters.  Applicable for any Lambda matrix, but slower than
##' other block-by-block methods.
##' Not exported.
##'
##' TODO:
##' (1) Write up quick note on theory (e.g. Laplace approximation).
##' (2) Figure out how to convert between full q-by-q matrix, and
##'     the format currently in the postVar attributes of the
##'     elements of the output of ranef.
##' (3) Test.
##' (4) Do we need to think carefully about the differences
##'     between REML and ML, beyond just multiplying by a different
##'     sigma^2 estimate?
##'
##' @param object \code{merMod} object
##' @return Sparse covariance matrix
condVar <- function(object) {
  s2 <- sigma(object)^2
  Lamt <- getME(object, "Lambdat")
  L <- getME(object, "L")

  ## never do it this way! fortune("SOOOO")
  #V <- solve(L, system = "A")
  #V <- chol2inv(L)
  #s2*crossprod(Lamt, V) %*% Lamt

  LL <- solve(L, Lamt, system = "A")
  s2 * crossprod(Lamt, LL)
}

mkMinimalData <- function(formula) {
    vars <- all.vars(formula)
    nVars <- length(vars)
    matr <- matrix(0, 2, nVars)
    data <- as.data.frame(matr)
    setNames(data, vars)
}

##' Make template for mixed model parameters
mkParsTemplate <- function(formula, data){
    if(missing(data)) data <- mkMinimalData(formula)
    mfRanef <- model.frame( subbars(formula), data)
    mmFixef <- model.matrix(nobars(formula) , data)
    reTrms <- mkReTrms(findbars(formula), mfRanef)
    cnms <- reTrms$cnms
    thetaNamesList <- mapply(mkPfun(), names(cnms), cnms)
    thetaNames <- unlist(thetaNamesList)
    betaNames <- colnames(mmFixef)
    list(beta  = setNames(numeric(length( betaNames)),  betaNames),
         theta = setNames(reTrms$theta, thetaNames),
         sigma = 1)
}

##' Make template for mixed model data
##'
##' Useful for simulating balanced designs and for
##' getting started on unbalanced simulations
##'
##' @param formula formula
##' @param data data -- not necessary
##' @param nGrps number of groups per grouping factor
##' @param rfunc function for generating covariate data
##' @param ... additional parameters for rfunc
mkDataTemplate <- function(formula, data,
                           nGrps = 2, nPerGrp = 1,
                           rfunc = NULL, ...){
    if(missing(data)) data <- mkMinimalData(formula)
    grpFacNames <- unique(barnames(findbars(formula)))
    varNames <- all.vars(formula)
    covariateNames <- setdiff(varNames, grpFacNames)
    nGrpFac <- length(grpFacNames)
    nCov <- length(covariateNames)
    grpFac <- gl(nGrps, nPerGrp)
    grpDat <- expand.grid(replicate(nGrpFac, grpFac, simplify = FALSE))
    colnames(grpDat) <- grpFacNames
    nObs <- nrow(grpDat)
    if(is.null(rfunc)) rfunc <- function(n, ...) rep(0, n)
    params <- c(list(nObs), list(...))
    covDat <- as.data.frame(replicate(nCov, do.call(rfunc, params),
                                      simplify = FALSE))
    colnames(covDat) <- covariateNames
    cbind(grpDat, covDat)
}

##' very flexible and convenient wrt formula,
##' very unflexible wrt everything else
##'
##' starting to get a little too sugary?
quickSimulate <- function(formula, nGrps, nPerGrp, family = gaussian) {
    pr <- mkParsTemplate(formula)
    dt <- mkDataTemplate(formula, nGrps = nGrps, nPerGrp = nPerGrp, rfunc = rnorm)
    response <- deparse(formula[[2]])
    dt[[response]] <- simulate(formula, newdata = dt, newparams = pr, family = family)[[1]]
    return(dt)
}

#----------------------------------------------------------------------
# formula parsing sugar
#----------------------------------------------------------------------

##' these functions pick up where findbars leaves off, in terms of sugar
##' @param REtrm an element of the result of findbars
##' @param REtrms the result of findbars
##' @return \code{reexpr} gives a one-sided formula with the linear
##' model formula for the raw model matrix. \code{grpfact} gives an
##' expression with the name of the grouping factor associated with
##' the raw model matrix. \code{termnms} gives a character vector with
##' the names of the random effects terms.
reexpr <- function(REtrm) substitute( ~ foo, list(foo = REtrm[[2]]))
grpfact <- function(REtrm) substitute(factor(fac), list(fac = REtrm[[3]]))
termnms <- function(REtrms) vapply(REtrms, safeDeparse, "")

##' mmList(): list of model matrices
##' ------    called from getME() & model.matrix(*, "randomListRaw")
mmList <- function(object, ...) UseMethod("mmList")
mmList.merMod <- function(object, ...) mmList(formula(object), model.frame(object))
mmList.formula <- function(object, frame, ...) {
    bars <- findbars(object)
    mm <- setNames(lapply(bars, function(b) model.matrix(eval(reexpr(b), frame), frame)),
		   termnms(bars))
    grp <- lapply(lapply(bars, grpfact), eval, frame)
    nl <- vapply(grp, nlevels, 1L)
    if (any(diff(nl) > 0))
        mm[order(nl, decreasing=TRUE)]
    else
        mm
}
##' examples  ---FIXME?--- put in tests // or export + 'real examples'
if(FALSE) {
library(lme4)
m <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
gm <- glmer(cbind(incidence, size-incidence) ~ period + (1|herd), cbpp, binomial)
simForm <- y ~ x + z + (x | f) + (z | g)
simDat <- lme4:::quickSimulate(simForm, 10, 5)
simDat <- simDat[simDat$f != "10", ] # unbalancedish design requiring
                                     # a flip in the order of terms
sm <- lmer(simForm, simDat)
lme4:::mmList.merMod(m)
lme4:::mmList.merMod(gm)
smmm <- lme4:::mmList.merMod(sm)
}

nloptwrap <- local({
    ## define default control values in environment of function ...
    defaultControl <- list(algorithm="NLOPT_LN_BOBYQA",
                           xtol_abs=1e-6, ftol_abs=1e-6, maxeval=1e5)
    ##
    function(par,fn,lower,upper,control=list(),...) {
        for (n in names(defaultControl))
            if (is.null(control[[n]])) control[[n]] <- defaultControl[[n]]
        res <- nloptr(x0=par, eval_f=fn, lb=lower,ub=upper, opts=control, ...)
        with(res,list(par=solution,
                      fval=objective,
                      feval=iterations,
                      conv=if (status>0) 0 else status,
                      message=message))
    }
})

nlminbwrap <- function(par, fn, lower, upper, control=list(), ...) {
    if (!is.null(control$maxfun)) {
        control$eval.max <- control$maxfun
        control$maxfun <- NULL
    }
    res <- nlminb(start = par, fn, gradient = NULL, hessian = NULL,
                  scale = 1, lower = lower, upper = upper,
                  control = control, ...)
    list(par = res$par, fval = res$objective,
         conv = res$convergence, message = res$message)
}

glmerLaplaceHandle <- function(pp, resp, nAGQ, tol, maxit, verbose) {
    .Call(glmerLaplace, pp, resp, nAGQ, tol, as.integer(maxit), verbose)
}

isFlexLambda <- function() FALSE