R/spde.common.R

Defines functions inla.dBind inla.extract.el inla.regex.match inla.extract.el.matrix inla.extract.el.data.frame inla.extract.el.list inla.spde.homogenise_B_matrix inla.matern.cov inla.matern.cov.s2 inla.spde.models inla.spde.sample inla.spde.sample.default inla.spde.sample.inla.spde inla.spde.precision inla.spde.result inla.spde.make.index inla.row.kron inla.spde.make.block.A inla.spde.make.A rbind.inla.data.stack.info inla.stack.remove.unused inla.stack.compress inla.stack inla.stack.sum inla.stack.join inla.stack.index inla.stack.do.extract inla.stack.LHS inla.stack.RHS inla.stack.data inla.stack.A

Documented in inla.dBind inla.extract.el inla.extract.el.data.frame inla.extract.el.list inla.extract.el.matrix inla.matern.cov inla.matern.cov.s2 inla.row.kron inla.spde.make.A inla.spde.make.block.A inla.spde.make.index inla.spde.models inla.spde.precision inla.spde.result inla.spde.sample inla.spde.sample.default inla.spde.sample.inla.spde inla.stack inla.stack.A inla.stack.compress inla.stack.data inla.stack.index inla.stack.join inla.stack.LHS inla.stack.remove.unused inla.stack.RHS inla.stack.sum rbind.inla.data.stack.info

## Export: inla.dBind inla.extract.el inla.extract.el!data.frame
## Export: inla.extract.el!list inla.extract.el!matrix inla.matern.cov
## Export: inla.matern.cov.s2 inla.row.kron
## Export: inla.spde.make.A inla.spde.make.block.A inla.spde.make.index
## Export: inla.spde.models inla.spde.precision inla.spde.result
## Export: inla.spde.sample inla.spde.sample!default
## Export: inla.spde.sample!inla.spde
## Internal: inla.spde.homogenise_B_matrix inla.regex.match

inla.dBind <- function(...)
{
    .Deprecated("Matrix::bdiag")
    return(bdiag(...))
}

inla.extract.el <- function(M, ...)
{
    if (is.null(M))
        return(NULL)
    UseMethod("inla.extract.el", M)
}

inla.regex.match =  function(x, match) {
    return(strsplit(x, match)[[1]][1]=="")
}

inla.extract.el.matrix <- function(M, match, by.row=TRUE, ...)
{
    if (by.row) {
        return(M[sapply(rownames(M), inla.regex.match, match=match),,drop=FALSE])
    } else {
        return(M[,sapply(colnames(M), inla.regex.match, match=match),drop=FALSE])
    }
}

inla.extract.el.data.frame <- function(M, match, by.row=TRUE, ...)
{
    if (by.row) {
        return(M[sapply(rownames(M), inla.regex.match, match=match),,drop=FALSE])
    } else {
        return(M[,sapply(colnames(M), inla.regex.match, match=match),drop=FALSE])
    }
}

inla.extract.el.list <- function(M, match, ...)
{
    return(M[sapply(names(M), inla.regex.match, match=match)])
}



inla.spde.homogenise_B_matrix <- function(B, n.spde, n.theta)
{
    if (!is.numeric(B))
        stop("B matrix must be numeric.")
    if (is.matrix(B)) {
        if ((nrow(B) != 1) && (nrow(B) != n.spde)) {
            stop(inla.paste(list("B matrix has",
                                 as.character(nrow(B)),
                                 "rows but should have 1 or",
                                 as.character(n.spde),
                                 sep=" ")))
        }
        if ((ncol(B) != 1) && (ncol(B) != 1+n.theta)) {
            stop(inla.paste(list("B matrix has",
                                 as.character(ncol(B)),
                                 "columns but should have 1 or",
                                 as.character(1+n.theta),
                                 sep=" ")))
        }
        if (ncol(B) == 1) {
            return(cbind(as.vector(B), matrix(0.0, n.spde, n.theta)))
        } else if (ncol(B) == 1+n.theta) {
            if (nrow(B) == 1) {
                return(matrix(as.vector(B), n.spde, 1+n.theta, byrow=TRUE))
            } else if (nrow(B) == n.spde) {
                return(B)
            }
        }
    } else { ## !is.matrix(B)
        if ((length(B) == 1) || (length(B) == n.spde)) {
            return(cbind(B, matrix(0.0, n.spde, n.theta)))
        } else if (length(B) == 1+n.theta) {
            return(matrix(B, n.spde, 1+n.theta, byrow=TRUE))
        } else {
            stop(inla.paste(list("Length of B vector is",
                                 as.character(length(B)),
                                 "but should be 1,",
                                 as.character(1+n.theta), "or",
                                 as.character(n.spde)),
                            sep=" "))
        }
    }
    stop(inla.paste(list("Unrecognised structure for B matrix"),
                    sep=" "))
}



inla.matern.cov <- function(nu,kappa,x,d=1,corr=FALSE, norm.corr=FALSE, theta, epsilon=1e-8)
{
    if (missing(theta)) { ## Ordinary Matern
        y = kappa*abs(x)
        if (corr) {
            ok = (y>=epsilon)
            if (nu<=0) {
                covariance = y*0
                covariance[!ok] = 1-y/epsilon
            } else {
                covariance = y*0
                covariance[ok] =
                    2^(1-nu)/gamma(nu) * (y[ok])^nu*besselK(y[ok], nu)
                if (any(!ok)) {
                    scale = 2^(1-nu)/gamma(nu)
                    ## corr = scale y^nu K_nu(y)
                    ##      = 1 - b y^(2 nu) + o(y^(2 nu), 0 < \nu < 1
                    ##      = 1 - b y^2 + o(y^2), 1 <= \nu
                    ## (1-corr(eps)/vari)/eps^(2nu) = b
                    if (nu < 1) {
                        exponent = 2*nu
                    } else {
                        exponent = 2
                    }
                    corr.eps <-
                        scale * epsilon^nu * besselK(epsilon, nu)
                    b = (1-corr.eps)/epsilon^exponent
                    covariance[!ok] = (1-b*y[!ok]^exponent)
                }
            }
            return(covariance)
        } else {
            ok = (y>=epsilon)
            covariance = y*0
            covariance[ok] =
                2^(1-nu)/gamma(nu+d/2)/(4*pi)^(d/2)/kappa^(2*nu)*
                    (y[ok])^nu*besselK(y[ok], nu)
            if (any(!ok)) {
                if (nu>0) { ## Regular Matern case
                    vari = gamma(nu)/gamma(nu+d/2)/(4*pi)^(d/2)/kappa^(2*nu)
                    scale = 2^(1-nu)/gamma(nu)
                    ## corr = scale y^nu K_nu(y)
                    ##      = 1 - b y^(2 nu) + o(y^(2 nu), 0 < \nu < 1
                    ##      = 1 - b y^2 + o(y^2), 1 <= \nu
                    ## (1-corr(eps)/vari)/eps^(2nu) = b
                    if (nu < 1) {
                        exponent = 2*nu
                    } else {
                        exponent = 2
                    }
                    corr.eps <-
                        scale * epsilon^nu * besselK(epsilon, nu)
                    b = (1-corr.eps)/epsilon^exponent
                    covariance[!ok] = vari*(1-b*y[!ok]^exponent)
                } else if (nu==0) { ## Limiting Matern case
                    g = 0.577215664901484 ## Euler's constant
                    covariance[!ok] =
                        2/gamma(d/2)/(4*pi)^(d/2)*
                            (-log(y[!ok]/2)-g)
                } else { ## (nu<0)
                    ## TODO: check this...
                    covariance[!ok] =
                        ((2^(1-nu)/gamma(nu+d/2)/(4*pi)^(d/2)/kappa^(2*nu)*
                          gamma(nu)*2^(nu-1))*(1-(y[!ok]/epsilon)) +
                         (2^(1-nu)/gamma(nu+d/2)/(4*pi)^(d/2)/kappa^(2*nu)*
                          epsilon^nu*besselK(epsilon, nu))*(y[!ok]/epsilon))
                }
            }
            return(covariance)
        }
    } else { ## Oscillating covariances
        y = abs(x)
        if (d>2L) {
            warning('Dimension > 2 not implemented for oscillating models.')
        }
        freq.max = 1000/max(y)
        freq.n = 10000
        w = seq(0,freq.max,length.out=freq.n)
        dw = w[2]-w[1]
        spec = 1/(2*pi)^d/(kappa^4+2*kappa^2*cos(pi*theta)*w^2+w^4)^((nu+d/2)/2)
       if (d==1L) {
            covariance = y*0+spec[1]*dw
        } else {
            covariance = y*0
        }
        for (k in 2:freq.n) {
            if (d==1L) {
                covariance = covariance+2*cos(y*w[k])*spec[k]*dw
            } else {
                covariance = covariance + w[k]*besselJ(y*w[k],0)*spec[k]*dw
            }
        }

        if (norm.corr) {
            noise.variance = 1/covariance[1]
        } else {
            noise.variance = 1
        }

        return(covariance*noise.variance)
    }
}


inla.matern.cov.s2 <- function(nu,kappa,x,norm.corr=FALSE,theta=0)
{
    stopifnot(inla.require("orthopolynom"))
    y = cos(abs(x))

    freq.max = 40L
    freq.n = freq.max+1L
    w = 0L:freq.max
    spec = 1/(kappa^4+2*kappa^2*cos(pi*theta)*w*(w+1)+w^2*(w+1)^2)^((nu+1)/2)
    leg = orthopolynom::legendre.polynomials(freq.max)
    covariance = y*0
    for (k in 1:freq.n) {
        covariance = (covariance + (2*w[k]+1)/(4*pi)*spec[k]*
                      orthopolynom::polynomial.values(leg[k],y)[[1]])
    }

    if (norm.corr) {
        noise.variance = 1/covariance[1]
    } else {
        noise.variance = 1
    }

    return(covariance*noise.variance)
}



inla.spde.models <- function(function.names=FALSE)
{
    types = c("spde1", "spde2")
    models = list()
    for (t in types) {
        models[[t]] =
            do.call(what=paste("inla.", t, ".models", sep=""),
                    args=list())
        if (function.names) {
            models[[t]] = paste("inla.", t, ".", models[[t]], sep="")
        }
    }

    if (function.names) {
        models = as.vector(do.call(c, models))
    }

    return(models)
}


inla.spde.sample <- function(...)
{
    warning("inla.spde.sample is deprecated.  Please use inla.qsample() in combination with inla.spde.precision() instead.")
    UseMethod("inla.spde.sample")
}

inla.spde.sample.default =
    function(precision, seed=NULL, ...)
{
    return(inla.finn(precision,
                     seed=(inla.ifelse(is.null(seed),
                                       0L,
                                       seed)))$sample)
}

inla.spde.sample.inla.spde =
    function(spde, seed=NULL, ...)
{
    precision = inla.spde.precision(spde, ...)
    return(inla.spde.sample(precision, seed=seed))
}



inla.spde.precision <- function(...)
{
    UseMethod("inla.spde.precision")
}

inla.spde.result <- function(...)
{
    inla.require.inherits(list(...)[[1]], "inla", "First parameter")
    inla.require.inherits(list(...)[[2]], "character", "Second parameter")
    UseMethod("inla.spde.result", list(...)[[3]])
}







inla.spde.make.index <- function(name, n.spde, n.group=1, n.repl=1, ...)
{
    if ("n.mesh" %in% names(list(...))) {
        warning("'n.mesh' is deprecated, please use 'n.spde' instead.")
        if (missing(n.spde) || is.null(n.spde))
            n.spde = list(...)$n.mesh
    }

    name.group = paste(name, ".group", sep="")
    name.repl = paste(name, ".repl", sep="")
    out = list()
    out[[name]]       = rep(rep(1:n.spde, times=n.group), times=n.repl)
    out[[name.group]] = rep(rep(1:n.group, each=n.spde), times=n.repl)
    out[[name.repl]]  = rep(1:n.repl, each=n.spde*n.group)
    return(out)
}





inla.row.kron <- function(M1, M2, repl=NULL, n.repl=NULL, weights=NULL) {
    M1 = inla.as.dgTMatrix(M1)
    M2 = inla.as.dgTMatrix(M2)
    n = nrow(M1)
    if (is.null(repl)) {
        repl = rep(1L, n)
    }
    if (is.null(n.repl)) {
        n.repl = max(repl)
    }
    if (is.null(weights)) {
        weights = rep(1, n)
    } else if (length(weights)==1L) {
        weights = rep(weights[1], n)
    }

    if (FALSE) {
    ## Slow version:
        print(system.time({
            M = (sparseMatrix(i=numeric(0), j=numeric(0), x=integer(0),
                              dims=c(n, ncol(M1)*ncol(M2))))
            for (k in seq_len(n)) {
                M[k,] = kronecker(M1[k,,drop=FALSE], M2[k,,drop=FALSE])
            }
            M = inla.as.dgTMatrix(M)
            weights.ii = weights[1L + M@i]
            M = (sparseMatrix(i=(1L + M@i),
                              j=(1L + M@j + ncol(M)*(repl[M@i+1L]-1L)),
                      x=weights.ii*M@x,
                              dims=c(n, n.repl*ncol(M))))
        }))
        M.slow = M
    }

    ## Fast version:
    ## TODO: Check robustness for all-zero rows.
    ## TODO: Maybe move big sparseMatrix call outside the loop.
    ## TODO: Automatically choose M1 or M2 for looping.

##    print(system.time({
    n1 = (as.vector(sparseMatrix(i=1L+M1@i, j=rep(1L, length(M1@i)),
                                 x=1L, dims=c(n, 1))))
    n2 = (as.vector(sparseMatrix(i=1L+M2@i, j=rep(1L, length(M2@i)),
                                 x=1L, dims=c(n, 1))))

    M = (sparseMatrix(i=integer(0), j=integer(0), x=numeric(0),
                      dims=c(n, ncol(M2)*ncol(M1)*n.repl)))
    n1 = n1[1L+M1@i]
    for (k in unique(n1)) {
        sub = which(n1==k)
        n.sub = length(sub)

        i.sub = 1L+M1@i[sub]
        j.sub = 1L+M1@j[sub]
        o1 = order(i.sub, j.sub)
        jj = rep(seq_len(k), times=n.sub/k)

        i.sub = i.sub[o1]
        j.sub = (sparseMatrix(i=i.sub,
                              j=jj,
                              x=j.sub[o1],
                              dims=c(n, k)))
        x.sub = (sparseMatrix(i=i.sub,
                              j=jj,
                              x=weights[i.sub]*M1@x[sub][o1],
                              dims=c(n, k)))
        sub2 = which(is.element(1L+M2@i, i.sub))

        if (length(sub2) > 0) {
            i = 1L+M2@i[sub2]
            ii = rep(i, times=k)
            repl.i = repl[ii]

            M = (M +
                 sparseMatrix(i=ii,
                              j=(1L+rep(M2@j[sub2], times=k)+
                                 ncol(M2)*(as.vector(j.sub[i,])-1L)+
                                 ncol(M2)*ncol(M1)*(repl.i-1L)),
                              x=(rep(M2@x[sub2], times=k)*
                                 as.vector(x.sub[i,])),
                              dims=c(n, ncol(M2)*ncol(M1)*n.repl)))
        }
    }
##}))

    ## For debugging:
    ##    print(max(abs(M-M.slow)))

    ##    o2 = order(n2[1L+M2@i], M2@i, M2@j)

    return(M)
}




## Add A-matrix rows belonging to the same "block" with optional
## weights and optional rescaling, by dividing row i by a_i, for |a_i|>0:
## B(i) = {j; block(i)==block(j) and sum_k |A_jk| > 0 }
## "count":   a_i = #{j in B(i)} }
## "weights": a_i = \sum_{j in B(i)} weights_j
## "sum":     a_i = \sum_{j in B(i)} \sum_k A_jk weights_j
##
## This function makes use of the feature of sparseMatrix to sum all
## values for multiple instances of the same (i,j) pair.
inla.spde.make.block.A =
    function(A,
             block,
             n.block = max(block),
             weights = NULL,
             rescale = c("none", "count", "weights", "sum"))
{
    A = inla.as.dgTMatrix(A)
    N = nrow(A)
    ## length(block) should be == N or 1
    if (length(block) == 1L) {
        block = rep(block, N)
    }
    if (is.null(weights)) {
        weights = rep(1, N)
    }

    rescale = match.arg(rescale)

    if (!(rescale == "none")) {
        if (rescale == "count") {
            ## Count the non-zero rows within each block
            sums = (sparseMatrix(i = block,
                                 j = rep(1L, N),
                                 x = (rowSums(abs(A)) > 0) * 1.0,
                                 dims = c(n.block, 1L)
                                 ))[block]
        } else if (rescale == "weights") {
            ## Sum the weights within each block
            sums = (sparseMatrix(i = block,
                                 j = rep(1L, N),
                                 x = (rowSums(abs(A)) > 0) * weights,
                                 dims = c(n.block, 1L)
                                 ))[block]
        } else { ## (rescale == "sum"){
            ## Sum the weighted values within each block
            sums = (sparseMatrix(i = block,
                                 j = rep(1L, N),
                                 x = rowSums(A) * weights,
                                 dims = c(n.block, 1L)
                                 ))[block]
        }
        ## Normalise:
        ok = (abs(sums) > 0)
        weights[ok] = weights[ok]/sums[ok]
    }

    return(inla.as.dgTMatrix(sparseMatrix(i = block[1L+A@i],
                                          j = 1L+A@j,
                                          x = A@x * weights[1L+A@i],
                                          dims = c(n.block, ncol(A))
                                          )))
}



inla.spde.make.A =
    function(mesh = NULL,
             loc = NULL,
             index = NULL,
             group = NULL,
             repl = 1L,
             n.spde = NULL,
             n.group = NULL,
             n.repl = NULL,
             group.mesh = NULL,
             weights = NULL,
             A.loc = NULL,
             A.group = NULL,
             group.index = NULL,
             block = NULL,
             n.block = NULL,
             block.rescale = c("none", "count", "weights", "sum"),
             ...)
## Deprecated/obsolete parameters: n.mesh, group.method
{
    ## A.loc can be specified instead of mesh+loc, optionally with
    ## index supplied.
    ## A.group can be specified instead of group and/or group.mesh,
    ## optionally with group.index supplied.

    if ("n.mesh" %in% names(list(...))) {
        warning("'n.mesh' is deprecated, use 'n.spde' instead.")
        n.spde = list(...)$n.mesh
    }
    if ("group.method" %in% names(list(...))) {
        group.method =
            match.arg(list(...)$group.method, c("nearest", "S0", "S1"))
        warning(paste("'group.method=", group.method,
                      "' is deprecated.  Specify 'degree=",
                      switch(group.method, nearest="0", S0="0", S1="1"),
                      "' in inla.mesh.1d() instead.", sep=""))
    }

    if (is.null(mesh)) {
        if (is.null(A.loc) && is.null(n.spde))
            stop("At least one of 'mesh', 'n.spde', and 'A.loc' must be specified.")
        if (!is.null(A.loc)) {
            n.spde = ncol(A.loc)
        }
    } else if (inherits(mesh, "inla.spde")) {
        spde <- mesh
        mesh <- spde$mesh
        inla.require.inherits(spde$mesh, c("inla.mesh", "inla.mesh.1d"), "'spde$mesh'")
        n.spde = spde$n.spde
    } else {
        inla.require.inherits(mesh, c("inla.mesh", "inla.mesh.1d"), "'mesh'")
        if (inherits(mesh, "inla.mesh.1d")) {
            n.spde = mesh$m
        } else {
            n.spde = mesh$n
        }
    }
    if (!is.null(group.mesh)) {
        inla.require.inherits(group.mesh, "inla.mesh.1d", "'mesh'")
    }

    n.group =
        ifelse(!is.null(n.group),
               n.group,
               ifelse(!is.null(A.group),
                      nrow(A.group),
                      ifelse(!is.null(group.mesh),
                             group.mesh$m,
                             max(1,
                                 ifelse(is.null(group),
                                        1,
                                        ifelse(length(group)==0,
                                               1,
                                               max(group)))))))
    n.repl =
        ifelse(!is.null(n.repl),
               n.repl,
               max(1, ifelse(is.null(repl),
                             1,
                             ifelse(length(repl)==0,
                                    1,
                                    max(repl)))))


    ## Handle loc and index input semantics:
    if (is.null(loc)) {
        if (is.null(A.loc)) {
            A.loc = Diagonal(n.spde, 1)
        }
    } else {
        if (is.null(mesh))
            stop("'loc' specified but 'mesh' is NULL.")

        ## Handle loc given as SpatialPoints or SpatialPointsDataFrame object
        if (inherits(loc, "SpatialPoints") ||
            inherits(loc, "SpatialPointsDataFrame")) {
          loc <- inla.spTransform(coordinates(loc),
                                  CRS(proj4string(loc)),
                                  mesh$crs,
                                  passthrough=TRUE)
        }

        A.loc = inla.mesh.project(mesh, loc=loc)$A
    }
    if (is.null(index)) {
        index = seq_len(nrow(A.loc))
    }
    ## Now 'index' points into the rows of 'A.loc'

    if (is.null(n.block)) {
        n.block = ifelse(is.null(block), length(index), max(block))
    }
    block.rescale = match.arg(block.rescale)

    ## Handle group semantics:
    ## TODO: FIXME!!! group, group.index, group.mesh, A.group, etc
    if (!is.null(A.group)) {
        if (!is.null(group) || !is.null(group.mesh)) {
            warning("'A.group' has been specified; ignoring non-NULL 'group' or 'group.mesh'.")
        }
    } else if (!is.null(group.mesh)) {
        if (is.null(group)) {
            group = rep(group.mesh$mid[1], length(index))
        }
    } else if (is.null(group)) {
        group = rep(1L, length(index))
    } else if (length(group) == 1) {
        group = rep(group, length(index))
    }
    if (is.null(group.index)) {
        group.index = seq_len(length(group))
    }
    ## Now 'group.index' points into the rows of 'A.group' or 'group'
    if (length(group.index) != length(index)) {
        stop(paste("length(group.index) != length(index): ",
                   length(group.index), " != ", length(index),
                   sep=""))
    }

    if (!is.null(group.mesh) && is.null(A.group)) {
        A.group = inla.mesh.1d.A(group.mesh, loc=group)
    }
    ## Now 'group.index' points into the rows of 'A.group' or 'group'

    ## Handle repl semantics:
    if (is.null(repl)) {
        repl = rep(1, length(index))
    } else if (length(repl) == 1) {
        repl = rep(repl, length(index))
    } else if (length(repl) != length(index)) {
        stop(paste("length(repl) != length(index): ",
                   length(repl), " != ", length(index),
                   sep=""))
    }

    if (length(index) > 0L) {
        A.loc = inla.as.dgTMatrix(A.loc[index,,drop=FALSE])

        if (length(A.loc@i) > 0L) {
            if (is.null(weights)) {
                weights = rep(1, length(index))
            } else if (length(weights)==1L) {
                weights = rep(weights[1], length(index))
            }
            if (!is.null(block)) {
                ## Leave the rescaling until the block phase,
                ## so that the proper rescaling can be determined.
                block.weights = weights
                weights = rep(1, length(index))
            }

            if (!is.null(A.group)) {
                A.group = inla.as.dgTMatrix(A.group[group.index,,drop=FALSE])
                A = (inla.row.kron(A.group, A.loc,
                                   repl=repl, n.repl=n.repl,
                                   weights=weights))
                ## More general version:
                ## A = inla.row.kron(A.repl,
                ##                   inla.row.kron(A.group, A.loc),
                ##                   weights=weights))
            } else {
                i = 1L+A.loc@i
                group.i = group[group.index[i]]
                repl.i = repl[i]
                weights.i = weights[i]
                A = (sparseMatrix(i=i,
                                  j=(1L+A.loc@j+
                                     n.spde*(group.i-1L)+
                                     n.spde*n.group*(repl.i-1L)),
                                  x=weights.i*A.loc@x,
                                  dims=(c(length(index),
                                          n.spde*n.group*n.repl))))
            }
            if (!is.null(block)) {
                A = (inla.spde.make.block.A(A=A,
                                            block=block,
                                            n.block=n.block,
                                            weights=block.weights,
                                            rescale=block.rescale))
            }
        } else {
            A = (sparseMatrix(i=integer(0),
                              j=integer(0),
                              x=numeric(0),
                              dims=c(n.block, n.spde*n.group*n.repl)))
        }
    } else {
        A = (sparseMatrix(i=integer(0),
                          j=integer(0),
                          x=numeric(0),
                          dims=c(0L, n.spde*n.group*n.repl)))
    }
    return(A)
}






#' Data stacking for advanced INLA models
#' 
#' Functions for combining data, effects and observation matrices into
#' \code{inla.stack} objects, and extracting information from such
#' objects.
#' 
#' @param data A list or code{data.frame} of named data vectors.
#' Scalars are expanded to match the number of rows in the A matrices, or any non-scalar data vectors.
#' An error is given if the input is inconsistent.
#' @param A A list of observation matrices. Scalars are expanded to diagonal matrices matching the effect vector lengths.
#' An error is given if the input is inconsistent or ambiguous.
#' @param effects A collection of effects/predictors.  Each list element corresponds
#' to an observation matrix, and must either be a single vector, a
#' list of vectors, or a \code{data.frame}. Single-element effect vectors are expanded
#'  to vectors matching the number of columns in the corresponding A matrix.
#'  An error is given if the input is inconsistent or ombiguous.
#' @param tag A string specifying a tag for later identification.
#' @param compress If \code{TRUE}, compress the model by removing duplicated rows of
#' effects, replacing the corresponding A-matrix columns with a single column containing the sum.
#' @param remove.unused If \code{TRUE}, compress the model by removing rows of effects
#' corresponding to all-zero columns in the \code{A} matrix (and removing
#' those columns).
#' @param stack A \code{inla.data.stack} object, created by a call to \code{inla.stack},
#' \code{inla.stack.sum}, or \code{inla.stack.join}.
#' @param ... For \code{inla.stack.join}, two or more data stacks of class
#' \code{inla.data.stack}, created by a call to \code{inla.stack},
#' \code{inla.stack.sum}, or \code{inla.stack.join}.
#' For \code{inla.stack.data}, a list of variables to be joined with
#' the data list.
#' 
#' @details
#' For models with a single effects collection, the outer list container for \code{A} and \code{effects} may be omitted.
#' 
#' Component size definitions:
#' \itemize{
#'    \item[\eqn{n_l}{n_l}] effect blocks
#'    \item[\eqn{n_k}{n_k}] effects
#'    \item[\eqn{n_i}{n_i}] data values
#'    \item[\eqn{n_{j,l}}{n_jl}] effect size for block \eqn{l}{l}
#'    \item[\eqn{n_j}{n_j}] \eqn{= \sum_{l=1}^{n_l} n_{j,l}}{sum_l n_jl} total effect size
#' }
#' 
#'  Input:
#'  \itemize{
#'    \item[\code{data}] \eqn{(y^1, \ldots, y^p)}{} \eqn{p}{p} vectors,
#'    each of length \eqn{n_i}{n_i}
#'    \item[\code{A}] \eqn{(A^1, \ldots, A^{n_l})}{} matrices of size
#'    \eqn{n_i \times n_{j,l}}{n_i by n_jl}
#'    \item[\code{effects}] \eqn{\left((x^{1,1},\ldots,x^{n_k,1}), \ldots,
#'      (x^{1,n_l},\ldots,x^{n_k,n_l})\right)}{} collections of effect vectors
#'    of length \eqn{n_{j,l}}{n_jl}
#'  }
#'  
#'  \deqn{
#'    \mbox{predictor}(y^1, \ldots, y^p) \sim
#'    \sum_{l=1}^{n_l} A^l \sum_{k=1}^{n_k} g(k, x^{k,l})
#'    = \tilde{A} \sum_{k=1}^{n_k} g(k, \tilde{x}^k)
#'  }{
#'    predictor(y^1, \ldots, y^p)
#'    ~ sum_{l=1}^{n_l} A^l sum_{k=1}^{n_k} g(k, x^{k,l})
#'    = tilde{A} sum_{k=1}^{n_k} g(k, tilde{x}^k)
#'  }
#'  where
#'  \deqn{
#'    \tilde{A} = \mbox{cbind}\left( A^1, \ldots, A^{n_l} \right)
#'  }{
#'    tilde{A} = cbind( A^1, ..., A^{n_l} )
#'  }
#'  \deqn{
#'    \tilde{x}^k = \mbox{rbind}\left( x^{k,1}, \ldots, x^{k,n_l} \right)
#'  }{
#'    tilde{x}^k = rbind( x^{k,1}, ..., x^{k,n_l} )
#'  }
#'  and for each block \eqn{l}{l}, any missing \eqn{x^{k,l}} is replaced by an
#'  \code{NA} vector.
#'
#' @return A data stack of class \code{inla.data.stack}.  Elements:
#'  \itemize{
#'    \item{\code{data} }{\eqn{=(y^1, \ldots, y^p, \tilde{x}^1, \ldots, \tilde{x}^{n_k})}}
#'    \item{\code{A} }{\eqn{=\tilde{A}}}
#'    \item{\code{data.names} }{List of data names, length \eqn{p}}
#'    \item{\code{effect.names} }{List of effect names, length \eqn{n_k}}
#'    \item{\code{n.data} }{Data length, \eqn{n_i}}
#'    \item{\code{index} }{List indexed by \code{tag}s, each element indexing
#'    into \eqn{i=1, \ldots, n_i}}
#'  }
#'  
#' @seealso 
#'  \code{\link{inla.spde.make.A}},
#'  \code{\link{inla.spde.make.index}}
#'
#' @keywords fmesher
#'  
#' @name inla.stack
#'
#' @examples
#' n <- 200
#' loc <- matrix(runif(n*2), n, 2)
#' mesh <- inla.mesh.2d(loc.domain = loc,
#'                      max.edge=c(0.05, 0.2))
#' proj.obs <- inla.mesh.projector(mesh, loc = loc)
#' proj.pred <- inla.mesh.projector(mesh, loc = mesh$loc)
#' spde <- inla.spde2.pcmatern(mesh,
#'                             prior.range = c(0.01, 0.01),
#'                             prior.sigma = c(10, 0.01))
#' 
#' covar <- rnorm(n)
#' field <- inla.qsample(n = 1, Q = inla.spde.precision(spde, theta = log(c(0.5, 1))))[,1]
#' y <- 2*covar + inla.mesh.project(proj.obs, field)
#' 
#' A.obs <- inla.spde.make.A(mesh, loc = loc)
#' A.pred = inla.spde.make.A(mesh, loc = proj.pred$loc)
#' stack.obs <-
#'     inla.stack(data=list(y=y),
#'                A=list(A.obs, 1),
#'                effects=list(c(list(Intercept = 1),
#'                               inla.spde.make.index("spatial", spde$n.spde)),
#'                             covar=covar),
#'                tag="obs")
#' stack.pred <-
#'     inla.stack(data=list(y=NA),
#'                A=list(A.pred),
#'                effects=list(c(list(Intercept = 1),
#'                               inla.spde.make.index("spatial", mesh$n))),
#'                tag="pred")
#' stack <- inla.stack(stack.obs, stack.pred)
#' 
#' formula <- y ~ -1 + Intercept + covar + f(spatial, model=spde)
#' result1 <- inla(formula,
#'                 data=inla.stack.data(stack.obs, spde = spde),
#'                 family="gaussian",
#'                 control.predictor = list(A = inla.stack.A(stack.obs),
#'                                         compute = TRUE))
#' 
#' plot(y, result1$summary.fitted.values[inla.stack.index(stack.obs,"obs")$data, "mean"],
#'      main = "Observations vs posterior predicted values at the data locations")
#' 
#' result2 <- inla(formula,
#'                 data=inla.stack.data(stack, spde = spde),
#'                 family="gaussian",
#'                 control.predictor = list(A = inla.stack.A(stack),
#'                                          compute = TRUE))
#' 
#' field.pred <- inla.mesh.project(proj.pred,
#'   result2$summary.fitted.values[inla.stack.index(stack,"pred")$data, "mean"])
#' field.pred.sd <- inla.mesh.project(proj.pred,
#'   result2$summary.fitted.values[inla.stack.index(stack,"pred")$data, "sd"])
#' 
#' plot(field, field.pred, main = "True vs predicted field")
#' abline(0, 1)
#' image(inla.mesh.project(mesh,
#'                         field = field,
#'                         dims = c(200,200)),
#'       main = "True field")
#' image(inla.mesh.project(mesh,
#'                         field = field.pred,
#'                         dims = c(200,200)),
#'       main = "Posterior field mean")
#' image(inla.mesh.project(mesh,
#'                         field = field.pred.sd,
#'                         dims = c(200,200)),
#'       main = "Prediction standard deviation")
#' plot(field, (field.pred - field) / 1,
#'      main = "True field vs standardised prediction residuals")
NULL



#' Internal function for merging raw stack information
#' 
#' @method rbind inla.data.stack.info
#' @export
rbind.inla.data.stack.info <- function(...)
{
    l = list(...)
    names(l) = NULL
    names.tmp = do.call(c, lapply(l, function(x) x$names))
    ncol.tmp = do.call(c, lapply(l, function(x) x$ncol))

    ncol = c()
    names = list()
    for (k in 1:length(names.tmp)) {
        name = names(names.tmp)[k]
        if (!is.null(names[[name]])) {
            if (!identical(names[[name]],
                           names.tmp[[k]])) {
                stop("Name mismatch.")
            }
        }
        names[[name]] = names.tmp[[k]]

        if (!is.null(as.list(ncol)[[name]])) {
            if (ncol[name] != ncol.tmp[[k]]) {
                stop("ncol mismatch.")
            }
        }
        ncol[name] = ncol.tmp[[k]]
    }

    external.names = names(names)
    internal.names = do.call(c, names)

    factors = rep(FALSE, length(internal.names))
    names(factors) = internal.names
    factor.names <-
        lapply(l, function(x) do.call(c,
                                      x$names[do.call(c,
                                                      lapply(x$data,
                                                             is.factor))]))
    for (factor.loop in seq_along(l)) {
        factors[factor.names[[factor.loop]]] = TRUE
    }

    handle.missing.columns <- function(x) {
        missing.names =
            setdiff(internal.names,
                    do.call(c, x$names))
        if (length(missing.names)>0) {
            df <- c(rep(list(rep(NA, x$nrow)),
                        sum(!factors[missing.names])),
                    rep(list(rep(as.factor(NA), x$nrow)),
                        sum(factors[missing.names])))
            names(df) <- c(missing.names[!factors[missing.names]],
                           missing.names[factors[missing.names]])
            df = as.data.frame(df)
            return(cbind(x$data, df))
        } else {
            return(x$data)
        }
    }

    data = do.call(rbind, lapply(l, handle.missing.columns))

    offset = 0
    index = list()
    for (k in 1:length(l)) {
        for (j in 1:length(l[[k]]$index)) {
            if (is.null(index[[names(l[[k]]$index)[j]]])) {
                index[[names(l[[k]]$index)[j]]] = l[[k]]$index[[j]] + offset
            } else {
                index[[names(l[[k]]$index)[j]]] =
                    c(index[[names(l[[k]]$index)[j]]],
                      l[[k]]$index[[j]] + offset)
            }
        }
        offset = offset + l[[k]]$nrow
    }

    info =
        list(data=data,
             nrow=nrow(data),
             ncol=ncol,
             names=names,
             index=index)
    class(info) = "inla.data.stack.info"

    return(info)
}

#' @describeIn inla.stack Remove unused entries from an existing stack
#' 
#' @export
inla.stack.remove.unused <- function(stack)
{
    inla.require.inherits(stack, "inla.data.stack", "'stack'")

    if (stack$effects$nrow<2) {
        return(stack)
    }

    ## Remove components with no effect:
    remove = rep(FALSE, stack$effects$nrow)
    remove.unused.indices =
        which(colSums(abs(stack$A[,,drop=FALSE]))==0)
    remove[remove.unused.indices] = TRUE

    index.new = rep(as.integer(NA), stack$effect$nrow)

    ncol.A = sum(!remove)
    if (ncol.A>0)
        index.new[!remove] = 1:ncol.A
    index.new[remove] = index.new[index.new[remove]]

    for (k in 1:length(stack$effects$index)) {
        stack$effects$index[[k]] = index.new[stack$effects$index[[k]]]
    }

    A = inla.as.dgTMatrix(stack$A)
    j.new = index.new[A@j+1L]
    ## Check for any zero-elements in remove.unused-columns:
    ok = !is.na(j.new)
    stack$A =
        sparseMatrix(i=A@i[ok]+1L,
                     j=j.new[ok],
                     x=A@x[ok],
                     dims=c(nrow(A), ncol.A))

    stack$effects$data = stack$effects$data[!remove,, drop=FALSE]
    stack$effects$nrow = ncol.A

    return(stack)
}

#' @describeIn inla.stack Compress an existing stack by removing duplicates
#'
#' @export
inla.stack.compress <- function(stack, remove.unused=TRUE)
{
    inla.require.inherits(stack, "inla.data.stack", "'stack'")

    if (stack$effects$nrow<2) {
        return(stack)
    }

    ii = do.call(order, as.list(stack$effects$data))
    jj.dupl =
        which(1L==
              diff(c(duplicated(stack$effects$data[ii,,drop=FALSE]),
                     FALSE)))
    kk.dupl =
        which(-1L==
              diff(c(duplicated(stack$effects$data[ii,,drop=FALSE]),
                     FALSE)))
    ## ii[jj.dupl] are the rows that have duplicates.
    ## ii[(jj.dupl[k]+1):kk.dupl[k]] are the duplicate rows for each k

    remove = rep(FALSE, stack$effects$nrow)
    index.new = rep(as.integer(NA), stack$effect$nrow)

    if (length(jj.dupl)>0) {
        for (k in 1:length(jj.dupl)) {
            i = ii[jj.dupl[k]]
            j = ii[(jj.dupl[k]+1):kk.dupl[k]]

            remove[j] = TRUE
            index.new[j] = i
        }
    }

    ncol.A = sum(!remove)
    if (ncol.A>0)
        index.new[!remove] = 1:ncol.A
    index.new[remove] = index.new[index.new[remove]]

    for (k in 1:length(stack$effects$index)) {
        stack$effects$index[[k]] = index.new[stack$effects$index[[k]]]
    }

    A = inla.as.dgTMatrix(stack$A)
    j.new = index.new[A@j+1L]
    ## Check for any zero-elements in remove.unused-columns:
    ok = !is.na(j.new)
    stack$A =
        sparseMatrix(i=A@i[ok]+1L,
                     j=j.new[ok],
                     x=A@x[ok],
                     dims=c(nrow(A), ncol.A))

    stack$effects$data = stack$effects$data[!remove,, drop=FALSE]
    stack$effects$nrow = ncol.A

    if (remove.unused) {
        return(inla.stack.remove.unused(stack))
    } else {
        return(stack)
    }
}



#' @describeIn inla.stack Shorthand for inla.stack.join and inla.stack.sum
#'
#' @export
inla.stack <- function(..., compress=TRUE, remove.unused=TRUE)
{
    if (all(sapply(list(...), function(x) inherits(x, "inla.data.stack")))) {
        return(do.call(inla.stack.join,
                       c(list(...),
                         compress=compress,
                         remove.unused=remove.unused)))
    } else {
        return(do.call(inla.stack.sum,
                       c(list(...),
                         compress=compress,
                         remove.unused=remove.unused)))
    }
}


#' @describeIn inla.stack Create data stack as a sum of predictors
#'
#' @export
inla.stack.sum <- function(data, A, effects,
                           tag="",
                           compress=TRUE,
                           remove.unused=TRUE)
{
    input.nrow <- function(x) {
        return(inla.ifelse(is.matrix(x) || is(x, "Matrix"),
                           nrow(x),
                           inla.ifelse(is.data.frame(x),
                                       rep(nrow(x), ncol(x)),
                                       length(x))))
    }
    input.ncol <- function(x) {
        return(inla.ifelse(is.matrix(x) || is(x, "Matrix"),
                           ncol(x),
                           inla.ifelse(is.data.frame(x),
                                       rep(1L, ncol(x)),
                                       1L)))
    }

    input.list.nrow <- function(l) {
        if (is.data.frame(l))
            return(input.nrow(l))
        return(do.call(c, lapply(l, input.nrow)))
    }
    input.list.ncol <- function(l) {
        if (is.data.frame(l))
            return(input.ncol(l))
        return(do.call(c, lapply(l, input.ncol)))
    }
    input.list.names <- function(l) {
        if (is.data.frame(l))
            return(colnames(l))
        is.df = sapply(l, is.data.frame)
        name = vector("list", length(l))
        if (!is.null(names(l)))
            name[!is.df] =
                lapply(names(l)[!is.df],
                       function(x) list(x))
        else
            name[!is.df] = ""
        name[is.df] =
            lapply(l[is.df],
                   function(x) as.list(colnames(x)))

        return(do.call(c, name))
    }


    parse.input.list <- function(l, n.A, error.tag, tag="", n.A.strict = FALSE) {
        ncol = input.list.ncol(l)
        nrow = input.list.nrow(l)
        names = input.list.names(l)
        if ((n.A>1) && any(nrow==1)) {
            for (k in which(nrow==1)) {
                if (ncol[k]==1) {
                    l[[k]] = rep(l[[k]], n.A)
                    nrow[k] = n.A
                } else {
                    stop(paste(error.tag,
                               "Automatic expansion only available for scalars.",
                               sep=""))
                }
            }
        }

        if (length(unique(c(names, ""))) < length(c(names, ""))) {
            stop(paste(error.tag,
                       "All variables must have unique names\n",
                       "Names: ('",
                       paste(names, collapse="', '", sep=""),
                       "')",
                       sep=""))
        }

        for (k in 1:length(names)) {
            if (ncol[k]==1) {
                names(names)[k] = names[[k]][[1]]
                names[[k]] = c(names[[k]][[1]])
            } else {
                names(names)[k] = names[[k]][[1]]
                names[[k]] = paste(names[[k]][[1]], ".", 1:ncol[k], sep="")
            }
        }

        names(nrow) = names(names)
        names(ncol) = names(names)

        ## data = as.data.frame(do.call(cbind, l))
        data = as.data.frame(l)
        names(data) = do.call(c, names)
        nrow = nrow(data)
        if ((n.A>1 || n.A.strict) && (nrow != n.A)) {
            stop(paste(error.tag,
                       "Mismatching row sizes: ",
                       paste(nrow, collapse=",", sep=""),
                       ", n.A=", n.A,
                       sep=""))
        }

        index = list(1:nrow)
        if (!is.null(tag)) {
            names(index) = tag
        }

        info = list(data=data, nrow=nrow, ncol=ncol, names=names, index=index)
        class(info) = "inla.data.stack.info"

        return(info)
    }

    if (is.null(tag))
        stop("'tag' must not be 'NULL'")

    ## Check if only a single block was specified.
    if (!is.list(A)) {
        A = list(A)
        effects = list(effects)
    }
    if (length(A) != length(effects))
        stop(paste("length(A)=", length(A),
                   " should be equal to length(effects)=", length(effects), sep=""))

    n.effects = length(effects)

    eff = list()
    for (k in 1:n.effects) {
        if (is.data.frame(effects[[k]])) {
            eff[[k]] =
                parse.input.list(list(effects[[k]]),
                                 input.ncol(A[[k]]),
                                 paste("Effect block ", k, ":\n", sep=""),
                                 tag)
        } else {
            if (!is.list(effects[[k]])) {
                tmp =
                    inla.ifelse(is.null(names(effects)[k]),
                                "",
                                names(effects)[k])
                effects[[k]] = list(effects[[k]])
                names(effects[[k]]) = tmp
            }
            eff[[k]] =
                parse.input.list(effects[[k]],
                                 input.ncol(A[[k]]),
                                 paste("Effect block ", k, ":\n", sep=""),
                                 tag)
        }
    }

    for (k in 1:n.effects) {
        if (is.vector(A[[k]])) {
            A[[k]] = Matrix(A[[k]], input.nrow(A[[k]]), 1)
        }
        if ((input.ncol(A[[k]])==1) && (eff[[k]]$nrow>1)) {
            if (input.nrow(A[[k]])!=1)
                stop(paste("ncol(A) does not match nrow(effect) for block ",
                           k, ": ",
                           input.ncol(A[[k]]), " != ", eff[[k]]$nrow, sep=""))
            A[[k]] = Diagonal(eff[[k]]$nrow, A[[k]][1,1])
        } else if (input.ncol(A[[k]]) != eff[[k]]$nrow) {
            stop(paste("ncol(A) does not match nrow(effect) for block ",
                       k, ": ",
                       input.ncol(A[[k]]), " != ", eff[[k]]$nrow, sep=""))
        }
    }
    if (length(unique(input.list.nrow(A)))>1) {
        stop(paste("Row count mismatch for A: ",
                   paste(input.list.nrow(A), collapse=",", sep=""),
                   sep=""))
    }
    A.nrow = nrow(A[[1]])
    A.ncol = input.list.ncol(A)

    data =
        parse.input.list(inla.ifelse(is.data.frame(data),
                                     as.list(data),
                                     data),
                         A.nrow,
                         paste("Data block:\n", sep=""),
                         tag,
                         n.A.strict = TRUE)

    effects = do.call(rbind.inla.data.stack.info, eff)

    A.matrix = do.call(cbind, A)
    A.nrow = nrow(A.matrix)
    A.ncol = ncol(A.matrix)

    if (length(unique(c(names(data$names), names(effects$names)))) <
        length(c(names(data$names), names(effects$names)))) {
        stop(paste("Names for data and effects must not coincide.\n",
                   "Data names:   ",
                   paste(names(data$names), collapse=", ", sep=""),
                   "\n",
                   "Effect names: ",
                   paste(names(effects$names), collapse=", ", sep=""),
                   sep=""))
    }

    stack = list(A=A.matrix, data=data, effects=effects)
    class(stack) = "inla.data.stack"

    if (compress) {
        return(inla.stack.compress(stack, remove.unused=remove.unused))
    } else if (remove.unused) {
        return(inla.stack.remove.unused(stack))
    } else {
        return(stack)
    }
}

#' @describeIn inla.stack Join two or more data stacks
#'
#' @export
inla.stack.join <- function(..., compress=TRUE, remove.unused=TRUE)
{
    S.input = list(...)

    data <- do.call(rbind.inla.data.stack.info,
                    lapply(S.input, function(x) x$data))
    effects <- do.call(rbind.inla.data.stack.info,
                       lapply(S.input, function(x) x$effects))
    ## The .bdiag form of bdiag takes a list as input.
    A <- .bdiag(lapply(S.input, function(x) x$A))

    S.output = list(A=A, data=data, effects=effects)
    class(S.output) = "inla.data.stack"

    if (length(unique(c(names(data$names), names(effects$names)))) <
        length(c(names(data$names), names(effects$names)))) {
        stop(paste("Names for data and effects must not coincide.\n",
                   "Data names:   ",
                   paste(names(data$names), collapse=", ", sep=""),
                   "\n",
                   "Effect names: ",
                   paste(names(effects$names), collapse=", ", sep=""),
                   sep=""))
    }

    if (compress) {
        return(inla.stack.compress(S.output, remove.unused=remove.unused))
    } else if (remove.unused) {
        return(inla.stack.remove.unused(S.output))
    } else {
        return(S.output)
    }
}





#' @describeIn inla.stack Extract tagged indices
#'
#' @export
inla.stack.index <- function(stack, tag)
{
    inla.require.inherits(stack, "inla.data.stack", "'stack'")

    if (is.null(tag)) {
        return(list(data=as.vector(do.call(c, stack$data$index)),
                    effects=(stack$data$nrow+
                             as.vector(do.call(c, stack$effects$index)))))
    } else {
        return(list(data=as.vector(do.call(c, stack$data$index[tag])),
                    effects=(stack$data$nrow+
                             as.vector(do.call(c, stack$effects$index[tag])))))
    }
}

## Extract an object list from an inla.stack internal structure
inla.stack.do.extract <- function(dat)
{
    inla.require.inherits(dat, "inla.data.stack.info", "'dat'")

    handle.entry <- function(x)
    {
        if (dat$ncol[[x]]>1) {
            return(matrix(do.call(c,
                                  dat$data[dat$names[[x]]]),
                          dat$nrow,
                          dat$ncol[[x]]))
        } else if (is.factor(dat$data[[dat$names[[x]]]])) {
            return(dat$data[[dat$names[[x]]]])
        }
        return(as.vector(dat$data[[dat$names[[x]]]]))
    }

    out = lapply(names(dat$names), handle.entry)
    names(out) = names(dat$names)

    return(out)
}


#' @describeIn inla.stack Extract data associated with the "left hand side" of the model
#' (e.g. the data itself, \code{Ntrials}, \code{link}, \code{E})
#'
#' @export
inla.stack.LHS <- function(stack)
{
    inla.require.inherits(stack, "inla.data.stack", "'stack'")

    return(inla.stack.do.extract(stack$data))
}

#' @describeIn inla.stack Extract data associated with the "right hand side" of the model
#' (all the covariates/predictors)
#'
#' @export
inla.stack.RHS <- function(stack)
{
    inla.require.inherits(stack, "inla.data.stack", "'stack'")

    return(inla.stack.do.extract(stack$effects))
}

#' @describeIn inla.stack Extract data for an inla call, and optionally join with other variables
#' 
#' @export
inla.stack.data <- function(stack, ...)
{
    inla.require.inherits(stack, "inla.data.stack", "'stack'")

    return(c(inla.stack.do.extract(stack$data),
             inla.stack.do.extract(stack$effects),
             list(...)))
}

#' @describeIn inla.stack Extract the "A matrix" for control.predictor
#' 
#' @export
inla.stack.A <- function(stack)
{
    inla.require.inherits(stack, "inla.data.stack", "'stack'")
    return(stack$A)
}
inbo/INLA documentation built on Dec. 6, 2019, 9:51 a.m.