R/matrices.R

In lava: Latent Variable Models

###{{{

`matrices` <-
    function(x,...) UseMethod("matrices")

###{{{ matrices.lvm

mat.lvm <- function(x,ii=index(x),...) {
    A <- ii$A ## Matrix with fixed parameters and ones where parameters are free
    M1 <- ii$M1 ## Index of free and _unique_ regression parameters
    P <- ii$P  ## Matrix with fixed variance parameters and ones where parameters are free
    P1 <- ii$P1 ## Index of free and _unique_ regression parameters

    constrain.par <- names(constrain(x))
    parval <- list()


    parBelongsTo <- list(mean=seq_len(ii$npar.mean),
                         reg=seq_len(ii$npar.reg)+ii$npar.mean,
                         cov=seq_len(ii$npar.var)+ii$npar.mean+ii$npar.reg,
                         epar=seq_len(ii$npar.ex)+with(ii,npar.reg+npar.var+npar.mean),
                         cpar=numeric())

    idxA <- which(M1==1)
    pidxA <- parBelongsTo$reg
    if (ii$npar.reg>0) {
        A[idxA] <- pidxA
        for (p in ii$parname) {
            idx <- which((x$par==p))
            newval <- A[idx[1]]
            attributes(newval)$reg.idx <- idx
            attributes(newval)$reg.tidx <- which(t(x$par==p))
            parval[[p]] <- newval
            if (length(idx)>1) {
                idxA <- c(idxA,idx[-1])
                pidxA <- c(pidxA,rep(A[idx[1]],length(idx)-1))
                A[idx] <- A[idx[1]]
            }
        } ## duplicate parameters
    }

    pars.var <- parBelongsTo$cov
    idxdiag <- (seq(ncol(P1))-1)*ncol(P1) + seq(ncol(P1))
    idxP <- idxdiag[which(P1[idxdiag]==1)]
    pidxP <- pars.var[seq_len(length(idxP))]
    P[idxP] <- pidxP

    pars.off.diag <- pars.var
    if (length(pidxP)>0) {
        pars.off.diag <- pars.off.diag[-seq_len(length(pidxP))]
    }

    counter <- 0
    if (length(pars.off.diag)>0 & ncol(P)>1)
        for (i in seq_len(ncol(P1)-1))
            for (j in seq(i+1,nrow(P1))) {
                if (ii$P1[j,i]!=0) {
                    counter <- counter+1
                    pos <- c(j+(i-1)*ncol(P1),
                             i+(j-1)*ncol(P1))
                    P[j,i] <- P[i,j] <- pars.off.diag[counter]
                    idxP <- c(idxP,pos); pidxP <- c(pidxP,P[j,i],P[i,j])
                }
            }

    if (length(ii$covparname)>0)
        for (p in ii$covparname) {
            idx <- which(x$covpar==p)
            isOffDiag <- !(idx[1]%in%idxdiag)
            if (!(p%in%ii$parname)) {
                parval[[p]] <- P[idx[1]]
            }
            attributes(parval[[p]])$cov.idx <- idx
            if (length(idx)>1+isOffDiag) {
                P[idx[-seq(1+isOffDiag)]] <- parval[[p]]
            }
            if (ii$npar.reg>0 && p%in%ii$parname) {
                parBelongsTo$reg <- c(parBelongsTo$reg,p)
                idx.reg <- which(x$par==p)
                P[idx] <- A[idx.reg[1]]
                atr <- attributes(parval[[p]])
                parval[[p]] <- A[idx.reg[1]]
                attributes(parval[[p]]) <- atr
                idxP <- c(idxP,idx)
                pidxP <- c(pidxP,rep(P[idx[1]],length(idx)))
            } else {
                idxP <- c(idxP,idx[-seq(1+isOffDiag)])
                pidxP <- c(pidxP,rep(P[idx[1]],length(idx)-1-isOffDiag))
            }
        } ## duplicate parameters

    idxM <- c()
    pidxM <- seq_len(ii$npar.mean)
    v <- NULL

    ## named <- sapply(x$mean, function(y) is.character(y) & !is.na(y))
    fixed <- sapply(x$mean, function(y) is.numeric(y) & !is.na(y))
    v <- rep(0,length(x$mean))
    names(v) <- colnames(P)
    if (ii$npar.mean>0) {
        idxM <- which(ii$v1==1)
        v[idxM] <- pidxM
    }
    if (any(fixed))
        v[fixed] <- unlist(x$mean[fixed])

    for (p in ii$mparname) {
        idx <- which(x$mean==p)
        if (!(p%in%c(ii$parname,ii$covparname))) {
            if (length(idx)>1) {
                pidxM <- c(pidxM,rep(v[idx[1]],length(idx)-1))
                idxM <- c(idxM,idx[-1])
            }
            parval[[p]] <- v[idx[1]]
            v[idx] <- parval[[p]]
        }
        attributes(parval[[p]])$m.idx <- idx
        if (p %in% ii$covparname & !(p %in% ii$parname)) {
            parBelongsTo$cov <- c(parBelongsTo$cov,p)
            idx.2 <- which(x$covpar==p)
            v[idx] <- P[idx.2[1]]
            pidxM <- c(pidxM,rep(P[idx.2[1]],length(idx)))
            idxM <- c(idxM,idx)
        }
        if (p %in% ii$parname) {
            parBelongsTo$reg <- c(parBelongsTo$reg,p)
            idx.2 <- which(x$par==p)
            v[idx] <- A[idx.2[1]]
            pidxM <- c(pidxM,rep(A[idx.2[1]],length(idx)))
            idxM <- c(idxM,idx)
        }
    }

    ## Ex-parameters
    idxE <- NULL
    pidxE <- parBelongsTo$epar
    ## named <- sapply(x$exfix, function(y) is.character(y) & !is.na(y))
    fixed <- sapply(x$exfix, function(y) is.numeric(y) & !is.na(y))
    epar <- rep(0,length(x$exfix))
    names(epar) <- names(x$expar)
    if (!(ii$npar.ex==0)) {
        idxE <- which(ii$e1==1)
        epar[idxE] <- pidxE
    }
    if (any(fixed))
        epar[fixed] <- unlist(x$exfix[fixed])
    for (p in ii$eparname) {
        idx <- which(x$exfix==p)
        if (!(p%in%c(ii$parname,ii$covparname,ii$mparname))) {
            if (length(idx)>1) {
                idxE <- c(idxE,idx[-1])
                pidxE <- c(pidxE,rep(epar[idx[1]],length(idx)-1))
            }
            parval[[p]] <- epar[idx[1]]
        }
        attributes(parval[[p]])$e.idx <- idx

        if (length(idx)>1)
            epar[idx[-1]] <- parval[[p]]
        if (p %in% setdiff(ii$covparname,c(ii$parname,ii$mparname))) {
            parBelongsTo$cov <- c(parBelongsTo$cov,p)
            idx.2 <- which(x$covpar==p)
            epar[idx] <- P[idx.2[1]]
            pidxE <- c(pidxE,rep(P[idx.2[1]],length(idx)))
            idxE <- c(idxE,idx)
        }
        if (p %in% setdiff(ii$parname,ii$mparname)) {
            parBelongsTo$reg <- c(parBelongsTo$reg,p)
            idx.2 <- which(x$par==p)
            epar[idx] <- A[idx.2[1]]
            pidxE <- c(pidxE,rep(A[idx.2[1]],length(idx)))
            idxE <- c(idxE,idx)
        }
        if (p %in% ii$mparname) {
            parBelongsTo$mean <- c(parBelongsTo$mean,p)
            idx.2 <- which(x$mean==p)
            epar[idx] <- v[idx.2[1]]
            pidxE <- c(pidxE,rep(v[idx.2[1]],length(idx)))
            idxE <- c(idxE,idx)
        }
    }
    ee <- cbind(idxE,pidxE); rownames(ee) <- names(x$expar)[ee[,1]]

    ## Constrained...
    constrain.par <- names(constrain(x))
    constrain.idx <- NULL

    if (length(constrain.par)>0) {
        constrain.idx <- list()
        for (p in constrain.par) {
            reg.tidx <- reg.idx <- cov.idx <- m.idx <- e.idx <- NULL
            myc <- constrain(x)[[p]]
            xargs <- manifest(x)[na.omit(match(attributes(myc)$args,manifest(x)))]
            if (length(xargs)>0) {
                parval[xargs] <- 0
            }
            if (p%in%ii$parname.all) {
                reg.idx <- which(x$par==p)
                reg.tidx <- which(t(x$par==p))
            }
            if (p%in%ii$covparname.all) {
                cov.idx <- which(x$covpar==p)
            }
            if (p%in%ii$mparname.all) {
                m.idx <- which(x$mean==p)
            }
            if (p%in%ii$eparname.all) {
                e.idx <- which(x$exfix==p)
            }
            constrain.idx[[p]] <- list(reg.idx=reg.idx,reg.tidx=reg.tidx,cov.idx=cov.idx,m.idx=m.idx,e.idx=e.idx)
        }
    }

    parBelongsTo <- lapply(parBelongsTo,function(x) sort(unique(x)))


    return(list(mean=cbind(idxM,pidxM),
                reg=cbind(idxA,pidxA),
                cov=cbind(idxP,pidxP),
                epar=ee,
                parval=parval,
                constrain.idx=constrain.idx,
                parBelongsTo=parBelongsTo))

}



matrices.lvm <- function(x,pars,meanpar=NULL,epars=NULL,data=NULL,...) {
    ii <- index(x)
    pp <- c(rep(NA,ii$npar.mean),pars,epars)
    ##v <- NULL
    v <- ii$v0
    if (!is.null(meanpar) && length(meanpar)>0) {
        pp[seq(ii$npar.mean)] <- meanpar
        v[ii$mean[,1]] <- meanpar[ii$mean[,2]]
    }

    A <- ii$A
    A[ii$reg[,1]] <- pp[ii$reg[,2]]
    P <- ii$P
    P[ii$cov[,1]] <- pp[ii$cov[,2]]
    e <- NULL

    if (length(x$expar)>0) {
        e <- rep(0,length(x$expar))
        fixed <- sapply(x$exfix, function(y) is.numeric(y) & !is.na(y))
        if (any(fixed))
            e[fixed] <- unlist(x$exfix[fixed])
        if (nrow(ii$epar)>0)
            e[ii$epar[,1]] <- pp[ii$epar[,2]]
        names(e) <- names(x$expar)
    }

    parval <- lapply(ii$parval,function(x) {
        res <- pp[x];
        attributes(res) <- attributes(x);
        res })
    ## Constrained...
    constrain.par <- names(constrain(x))
    constrain.idx <- NULL
    cname <- constrainpar <- c()
    if (length(constrain.par)>0 && is.numeric(c(pars,meanpar))) {
        constrain.idx <- list()
        for (p in constrain.par) {
            cname <- c(cname,p)
            myc <- constrain(x)[[p]]
            xargs <- manifest(x)[na.omit(match(attributes(myc)$args,manifest(x)))]
            if (length(xargs)>0) {
                if (!is.null(data)) {
                    parval[xargs] <- (data)[xargs]
                } else parval[xargs] <- 0
            }
            val <- unlist(c(parval,constrainpar,x$mean,e)[attributes(myc)$args])
            cpar <- myc(val);
            constrainpar <- c(constrainpar,list(cpar)); names(constrainpar) <- cname
            if (p%in%ii$parname.all) {
                if (!is.null(val))
                    A[ii$constrain.idx[[p]]$reg.idx] <- cpar
            }
            if (p%in%ii$covparname.all) {
                if (!is.null(val))
                    P[ii$constrain.idx[[p]]$cov.idx] <- cpar
            }
            if (p%in%ii$mparname.all) {
                if (!is.null(val))
                    v[ii$constrain.idx[[p]]$m.idx] <- cpar
            }
            if (p%in%ii$eparname.all) {
                if (!is.null(val))
                    e[ii$constrain.idx[[p]]$e.idx] <- cpar
            }
        }
    }

    return(list(A=A, P=P, v=v, e=e, parval=parval,
                constrain.idx=ii$constrain.idx, constrainpar=constrainpar))
}

###}}} matrices.lvm

###{{{ matrices.multigroup

matrices.multigroup <- function(x, p, ...) {
    pp <- modelPar(x,p)
    res <- list()
    for (i in seq_len(x$ngroup))
        res <- c(res, list(matrices2(x$lvm[[i]],pp$p[[i]])))
    return(res)
}

###}}}

matrices2 <- function(x,p,...) {
    m0 <- p[seq_len(index(x)$npar.mean)]
    p0 <- p[with(index(x),seq_len(npar)+npar.mean)]
    e0 <- p[with(index(x),seq_len(npar.ex)+npar.mean+npar)]
    matrices(x,p0,m0,e0,...)
}

###{{{ matrices, to be superseeded by above definition

matrices.lvm <- function(x,pars,meanpar=NULL,epars=NULL,data=NULL,...) {
    ii <- index(x)
    A <- ii$A ## Matrix with fixed parameters and ones where parameters are free
    M1 <- ii$M1 ## Index of free and _unique_ regression parameters
    P <- ii$P  ## Matrix with fixed variance parameters and ones where parameters are free
    P1 <- ii$P1 ## Index of free and _unique_ regression parameters

    constrain.par <- names(constrain(x))
    parval <- list()

    if (ii$npar.reg>0) {
        A[which(M1==1)] <- pars[seq_len(ii$npar.reg)]
        for (p in ii$parname) {
            idx <- which((x$par==p))
            newval <- A[idx[1]]
            attributes(newval)$reg.idx <- idx
            attributes(newval)$reg.tidx <- which(t(x$par==p))
            parval[[p]] <- newval
            if (length(idx)>1) {
                A[idx[-1]] <- parval[[p]]
            }
        } ## duplicate parameters
    }

    if (ii$npar.reg==0) {
        pars.var <- pars
    } else {
        pars.var <- pars[-seq_len(ii$npar.reg)]
    }

    diag(P)[ii$which.diag] <- pars.var[seq_along(ii$which.diag)]
    pars.off.diag <- pars.var
    if (length(ii$which.diag)>0)
        pars.off.diag <- pars.off.diag[-seq_along(ii$which.diag)]
    counter <- 0
    if (length(pars.off.diag)>0 & ncol(P)>1)
        for (i in seq_len(ncol(P1)-1))
            for (j in seq(i+1,nrow(P1))) {
                if (ii$P1[j,i]!=0) {
                    counter <- counter+1
                    P[j,i] <- pars.off.diag[counter]
                }
            }

    if (length(ii$covparname)>0)
        for (p in ii$covparname) {
            idx <- which(x$covpar==p)
            if (!(p%in%ii$parname)) {
                parval[[p]] <- P[idx[1]]
            }
            attributes(parval[[p]])$cov.idx <- idx
            if (length(idx)>1) {
                P[idx[-1]] <- parval[[p]]
            }
            if (ii$npar.reg>0 && p%in%ii$parname) {
                idx.reg <- which(x$par==p)
                P[idx] <- A[idx.reg[1]]
                atr <- attributes(parval[[p]])
                parval[[p]] <- A[idx.reg[1]] ###?????
                attributes(parval[[p]]) <- atr
            }
        } ## duplicate parameters
    P[upper.tri(P)] <- t(P)[upper.tri(P)]  ## Symmetrize...


    v <- NULL
    {
        ## named <- sapply(x$mean, function(y) is.character(y) & !is.na(y))
        fixed <- sapply(x$mean, function(y) is.numeric(y) & !is.na(y))
        v <- rep(0,length(x$mean))
        names(v) <- colnames(P)
        if (!(is.null(meanpar) | ii$npar.mean==0))
            v[ii$v1==1] <- meanpar
        if (any(fixed))
            v[fixed] <- unlist(x$mean[fixed])

        for (p in ii$mparname) {
            idx <- which(x$mean==p)

            if (!(p%in%c(ii$parname,ii$covparname))) {
                parval[[p]] <- v[idx[1]]
            }
            attributes(parval[[p]])$m.idx <- idx

            if (length(idx)>1)
                v[idx[-1]] <- parval[[p]]
            if (p %in% ii$covparname & !(p %in% ii$parname)) {
                idx.2 <- which(x$covpar==p)
                v[idx] <- P[idx.2[1]]
            }
            if (p %in% ii$parname) {
                idx.2 <- which(x$par==p)
                v[idx] <- A[idx.2[1]]
            }
        }
    }


    ## Ex-parameters
    e <- NULL
    {
        ## named <- sapply(x$exfix, function(y) is.character(y) & !is.na(y))
        fixed <- sapply(x$exfix, function(y) is.numeric(y) & !is.na(y))

        e <- rep(0,length(x$exfix))

        names(e) <- names(x$expar)
        if (!(is.null(epars) | ii$npar.ex==0))
            e[which(ii$e1==1)] <- epars
        if (any(fixed))
            e[fixed] <- unlist(x$exfix[fixed])
        for (p in ii$eparname) {
            idx <- which(x$exfix==p)
            if (!(p%in%c(ii$parname,ii$covparname,ii$mparname))) {
                parval[[p]] <- e[idx[1]]
            }
            attributes(parval[[p]])$e.idx <- idx

            if (length(idx)>1)
                e[idx[-1]] <- parval[[p]]
            if (p %in% setdiff(ii$covparname,c(ii$parname,ii$mparname))) {
                idx.2 <- which(x$covpar==p)
                e[idx] <- P[idx.2[1]]
            }
            if (p %in% setdiff(ii$parname,ii$mparname)) {
                idx.2 <- which(x$par==p)
                e[idx] <- A[idx.2[1]]
            }
            if (p %in% ii$mparname) {
                idx.2 <- which(x$mean==p)
                e[idx] <- v[idx.2[1]]
            }
        }
    }

    ## Constrained...
    constrain.idx <- NULL
    cname <- constrainpar <- c()
    if (length(constrain.par)>0 && is.numeric(c(pars,meanpar,e))) {
        constrain.idx <- list()
        for (p in constrain.par) {
            cname <- c(cname,p)
            reg.tidx <- reg.idx <- cov.idx <- m.idx <- e.idx <- NULL
            myc <- constrain(x)[[p]]
            xargs <- manifest(x)[na.omit(match(attributes(myc)$args,manifest(x)))]
            if (length(xargs)>0) {
                if (!is.null(data)) {
                    parval[xargs] <- (data)[xargs]
                } else parval[xargs] <- 0
            }
            val <- rbind(unlist(c(parval,constrainpar,x$mean,e)[attributes(myc)$args]))
            cpar <- myc(val);
            constrainpar <- c(constrainpar,list(cpar)); names(constrainpar) <- cname
            if (p%in%ii$parname.all) {
                reg.idx <- which(x$par==p)
                reg.tidx <- which(t(x$par==p))
                if (!is.null(val))
                    A[reg.idx] <- cpar##myc(val)
            }
            if (p%in%ii$covparname.all) {
                cov.idx <- which(x$covpar==p)
                if (!is.null(val))
                    P[cov.idx] <- cpar##myc(val)
            }
            if (p%in%ii$mparname.all) {
                m.idx <- which(x$mean==p)
                if (!is.null(val))
                    v[m.idx] <- cpar##myc(val)
            }
            if (p%in%ii$eparname.all) {
                e.idx <- which(x$exfix==p)
                if (!is.null(val))
                    e[e.idx] <- cpar##myc(val)
            }
            constrain.idx[[p]] <- list(reg.idx=reg.idx,reg.tidx=reg.tidx,cov.idx=cov.idx,m.idx=m.idx,e.idx=e.idx)
        }
    }

    if (x$index$sparse & !is.character(class(pars)[1])) {
        A <- as(A,"sparseMatrix")
        P <- as(P,"sparseMatrix")
        v <- as(v,"sparseMatrix")
    }
    return(list(A=A, P=P, v=v, e=e, parval=parval, constrain.idx=constrain.idx, constrainpar=constrainpar))
}

###}}} matrices Obsolete