psindex: Parameter Stability Index

standardize.est.lv.x <- function(x, lavobject, partable = NULL, cov.std = TRUE) {
    # embed x in est
  
  
    est <- lav_object_inspect_est(lavobject)
    free.idx <- which(lavobject@ParTable$free > 0L)
    stopifnot(length(x) == length(free.idx))
    est[free.idx] <- x

    # take care of setResidualElements...
    lavmodel <- lav_model_set_parameters(lavmodel = lavobject@Model, x = x)
    GLIST <- lavmodel@GLIST

    standardize.est.lv(lavobject = lavobject, partable = partable, est = est, 
                       GLIST = GLIST, cov.std = cov.std)
}

standardize.est.all.x <- function(x, lavobject, partable = NULL, cov.std = TRUE) {
    # embed x in est
    est <- lav_object_inspect_est(lavobject)
    free.idx <- which(lavobject@ParTable$free > 0L)
    stopifnot(length(x) == length(free.idx))
    est[free.idx] <- x

    # take care of setResidualElements...
    lavmodel <- lav_model_set_parameters(lavmodel = lavobject@Model, x = x)
    GLIST <- lavmodel@GLIST

    standardize.est.all(lavobject = lavobject, partable = partable, est = est,
                        est.std = NULL, GLIST = GLIST, cov.std = cov.std)
}

standardize.est.all.nox.x <- function(x, lavobject, partable = NULL, cov.std = TRUE) {
    # embed x in est
    est <- lav_object_inspect_est(lavobject)
    free.idx <- which(lavobject@ParTable$free > 0L)
    stopifnot(length(x) == length(free.idx))
    est[free.idx] <- x

    # take care of setResidualElements...
    lavmodel <- lav_model_set_parameters(lavmodel = lavobject@Model, x = x)
    GLIST <- lavmodel@GLIST

    standardize.est.all.nox(lavobject = lavobject, partable = partable, est = est, 
                            est.std = NULL, GLIST = GLIST, cov.std = cov.std)
}

unstandardize.est.ov.x <- function(x, lavobject) {
    partable <- lavobject@ParTable
    partable$ustart <- x
    unstandardize.est.ov(partable=partable, ov.var=lavobject@SampleStats@var, 
                         cov.std=TRUE)
}

standardize.est.lv <- function(lavobject, partable=NULL, est=NULL, GLIST=NULL,
                               cov.std = TRUE) {

    if(is.null(partable)) partable <- lavobject@ParTable
    if(is.null(est))   est <- lav_object_inspect_est(lavobject)
    if(is.null(GLIST)) GLIST <- lavobject@Model@GLIST

    out <- est; N <- length(est)
    stopifnot(N == length(partable$lhs))

    nmat <- lavobject@Model@nmat

    # compute ETA
    LV.ETA <- computeVETA(lavmodel       = lavobject@Model,
                          GLIST          = GLIST)
    
    for(g in 1:lavobject@Model@nblocks) {

        ov.names <- vnames(lavobject@ParTable, "ov", block=g) # not user, 
                                                       # which may be incomplete
        lv.names <- vnames(lavobject@ParTable, "lv", block=g)
       
        # shortcut: no latents in this block, nothing to do
        if(length(lv.names) == 0L)
            next

        # which mm belong to block g?
        mm.in.group <- 1:nmat[g] + cumsum(c(0,nmat))[g]
        MLIST <- GLIST[ mm.in.group ]

        ETA2 <- diag(LV.ETA[[g]])
        ETA  <- sqrt(ETA2)

        # 1a. "=~" regular indicators
        idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) & 
                     partable$block == g)
        out[idx] <- out[idx] * ETA[ match(partable$lhs[idx], lv.names) ]

        # 1b. "=~" regular higher-order lv indicators
        idx <- which(partable$op == "=~" & !(partable$rhs %in% ov.names) &
                     partable$block == g)
        out[idx] <- ( out[idx] * ETA[ match(partable$lhs[idx], lv.names) ]
                               / ETA[ match(partable$rhs[idx], lv.names) ] )

        # 1c. "=~" indicators that are both in ov and lv
        #idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
        #                             & partable$rhs %in% lv.names &
        #             partable$block == g)

        # 2. "~" regressions (and "<~")
        idx <- which((partable$op == "~" | partable$op == "<~") & 
                     partable$lhs %in% lv.names &
                     partable$block == g)
        out[idx] <- out[idx] / ETA[ match(partable$lhs[idx], lv.names) ] 

        idx <- which((partable$op == "~" | partable$op == "<~") & 
                     partable$rhs %in% lv.names &
                     partable$block == g)
        out[idx] <- out[idx] * ETA[ match(partable$rhs[idx], lv.names) ]

        # 3a. "~~" ov
        #idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) & 
        #             partable$block == g)

        # 3b. "~~" lv
        # ATTENTION: in Mplus 4.1, the off-diagonal residual covariances 
        #            were computed by the formula cov(i,j) / sqrt(i.var*j.var)
        #            were i.var and j.var where diagonal elements of ETA
        #
        #            in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
        #            elements are the 'PSI' diagonal elements!!

        # variances
        rv.idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
                        partable$lhs == partable$rhs &
                        partable$block == g)
        out[rv.idx] <- ( out[rv.idx] / ETA[ match(partable$lhs[rv.idx], lv.names) ]
                                     / ETA[ match(partable$rhs[rv.idx], lv.names) ] )

        # covariances lv
        # three types:
        # - only lhs is LV (and fixed.x = FALSE)
        # - only rhs is LV (and fixed.x = FALSE)
        # - both lhs and rhs are LV (regular case)
        if(cov.std) {
            if(!is.complex(est[rv.idx])) {
                RV <- sqrt(abs(est[rv.idx])) # abs in case of heywood cases
            } else {
                RV <- sqrt( est[rv.idx] )
            }
            rv.names <- partable$lhs[rv.idx]
        }

        # left
        idx.lhs <- which(partable$op == "~~" &
                         partable$lhs %in% lv.names &
                         partable$lhs != partable$rhs &
                         partable$block == g)
        if(length(idx.lhs) > 0L) {
            if(cov.std == FALSE) {
                out[idx.lhs] <- 
                   (out[idx.lhs] / ETA[ match(partable$lhs[idx.lhs], lv.names)])
            } else {
                out[idx.lhs] <- 
                   (out[idx.lhs] / RV[ match(partable$lhs[idx.lhs], rv.names)])
            }
        }

        # right
        idx.rhs <- which(partable$op == "~~" & 
                         partable$rhs %in% lv.names &
                         partable$lhs != partable$rhs &
                         partable$block == g)
        if(length(idx.rhs) > 0L) {
            if(cov.std == FALSE) {
                out[idx.rhs] <- 
                    (out[idx.rhs] / ETA[ match(partable$rhs[idx.rhs],lv.names)])
            } else {
                out[idx.rhs] <- 
                    (out[idx.rhs] / RV[ match(partable$rhs[idx.rhs], rv.names)])
            }
        }


        # 4a. "~1" ov
        #idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
        #             partable$block == g)

        # 4b. "~1" lv
        idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
                     partable$block == g)
        out[idx] <- out[idx] / ETA[ match(partable$lhs[idx], lv.names) ]
    }

    # 5a ":="
    idx <- which(partable$op == ":=")
    if(length(idx) > 0L) {
        x <- out[ partable$free & !duplicated(partable$free) ]
        out[idx] <- lavobject@Model@def.function(x)
    }

    # 5b "=="
    idx <- which(partable$op == "==")
    if(length(idx) > 0L) {
        x <- out[ partable$free & !duplicated(partable$free) ]
        out[idx] <- lavobject@Model@ceq.function(x)
    }

    # 5c. "<" or ">"
    idx <- which((partable$op == "<" | partable$op == ">"))
    if(length(idx) > 0L) {
        x <- out[ partable$free & !duplicated(partable$free) ]
        out[idx] <- lavobject@Model@cin.function(x)
    }

    out
}

standardize.est.all <- function(lavobject, partable=NULL, est=NULL, est.std=NULL,
                                GLIST = NULL, cov.std = TRUE) {
# browser()
    if(is.null(partable)) partable <- lavobject@ParTable
    if(is.null(est)) est <- lav_object_inspect_est(lavobject)
    if(is.null(est.std)) {
        est.std <- standardize.est.lv(lavobject, partable = partable, est = est,
                                      GLIST = GLIST, cov.std = cov.std)
    }
    if(is.null(GLIST)) GLIST <- lavobject@Model@GLIST

    out <- est.std; N <- length(est.std)
    stopifnot(N == length(partable$lhs))

    VY <- computeVY(lavmodel = lavobject@Model, GLIST = GLIST,
                    diagonal.only = TRUE)


    for(g in 1:lavobject@Model@nblocks) {

        ov.names <- vnames(lavobject@ParTable, "ov", block = g) # not user
        lv.names <- vnames(lavobject@ParTable, "lv", block = g)

        OV2 <- VY[[g]]
        # replace zero values by NA (but keep negative values)
        zero.idx <- which(abs(OV2) < .Machine$double.eps)
        if(length(zero.idx) > 0L) {
            OV2[zero.idx] <- as.numeric(NA)
        }
     
        # replace negative values by NA (for sqrt)
        tmp.OV2 <- OV2
        neg.idx <- which(tmp.OV2 < 0)
        if(length(neg.idx) > 0L) {
            tmp.OV2[neg.idx] <- as.numeric(NA)
        }
        OV  <- sqrt(tmp.OV2)

        if(lavobject@Model@conditional.x) {
            # extend OV with ov.names.x
            ov.names.x <- vnames(lavobject@ParTable, "ov.x", block = g)
            ov.names.nox <- vnames(lavobject@ParTable, "ov.nox", block = g)
            ov.names <- c(ov.names.nox, ov.names.x)
            OV2 <- c(OV2, diag(lavobject@SampleStats@cov.x[[g]]))
            OV <- c(OV, sqrt(diag(lavobject@SampleStats@cov.x[[g]])))
        }

        # 1a. "=~" regular indicators
        idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
                     partable$block == g)
        out[idx] <- out[idx] / OV[ match(partable$rhs[idx], ov.names) ]

        # 1b. "=~" regular higher-order lv indicators

        # 1c. "=~" indicators that are both in ov and lv
        #idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
        #                             & partable$rhs %in% lv.names &
        #             partable$block == g)

        # 2. "~" regressions (and "<~")
        idx <- which((partable$op == "~" | partable$op == "<~") & 
                     partable$lhs %in% ov.names &
                     partable$block == g)
        out[idx] <- out[idx] / OV[ match(partable$lhs[idx], ov.names) ]

        idx <- which((partable$op == "~" | partable$op == "<~") & 
                     partable$rhs %in% ov.names &
                     partable$block == g)
        out[idx] <- out[idx] * OV[ match(partable$rhs[idx], ov.names) ]

        # 3a. "~~" ov
        # ATTENTION: in Mplus 4.1, the off-diagonal residual covariances 
        #            were computed by the formula cov(i,j) / sqrt(i.var*j.var)
        #            were i.var and j.var where diagonal elements of OV
        #
        #            in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
        #            elements are the 'THETA' diagonal elements!!

        # variances
        rv.idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) & 
                        partable$lhs == partable$rhs &
                        partable$block == g)
        #out[rv.idx] <- ( out[rv.idx] / OV[ match(partable$lhs[rv.idx], ov.names) ]
        #                             / OV[ match(partable$rhs[rv.idx], ov.names) ] )
        out[rv.idx] <- ( out[rv.idx] / 
                             OV2[ match(partable$lhs[rv.idx], ov.names) ] )

        # covariances ov
        # three types:
        # - only lhs is OV (and fixed.x = FALSE)
        # - only rhs is OV (and fixed.x = FALSE)
        # - both lhs and rhs are OV (regular case)
        if(cov.std) {
            if(!is.complex(est[rv.idx])) {
                RV <- sqrt(abs(est[rv.idx]))
            } else {
                RV <- sqrt( est[rv.idx] )
            }
            rv.names <- partable$lhs[rv.idx]
        }

        # left
        idx.lhs <- which(partable$op == "~~" &
                         !(partable$lhs %in% lv.names) &
                         partable$lhs != partable$rhs &
                         partable$block == g)
        if(length(idx.lhs) > 0L) {
            if(cov.std == FALSE) {
                out[idx.lhs] <- 
                   (out[idx.lhs] / OV[ match(partable$lhs[idx.lhs], ov.names)])
            } else {
                out[idx.lhs] <- 
                   (out[idx.lhs] / RV[ match(partable$lhs[idx.lhs], rv.names)])
            }
        }

        # right
        idx.rhs <- which(partable$op == "~~" & 
                         !(partable$rhs %in% lv.names) &
                         partable$lhs != partable$rhs &
                         partable$block == g)
        if(length(idx.rhs) > 0L) {
            if(cov.std == FALSE) {
                out[idx.rhs] <- 
                    (out[idx.rhs] / OV[ match(partable$rhs[idx.rhs], ov.names)])
            } else {
                out[idx.rhs] <- 
                    (out[idx.rhs] / RV[ match(partable$rhs[idx.rhs], rv.names)])
            }
        }

        # 3b. "~~" lv
        #idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
        #             partable$block == g)

        # 4a. "~1" ov
        idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
                     partable$block == g)
        out[idx] <- out[idx] / OV[ match(partable$lhs[idx], ov.names) ]

        # 4b. "~1" lv
        #idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
        #             partable$block == g)

        # 4c. "|" thresholds
        idx <- which(partable$op == "|" & !(partable$lhs %in% lv.names) &
                     partable$block == g)
        out[idx] <- out[idx] / OV[ match(partable$lhs[idx], ov.names) ]

        # 4d. "~*~" scales
        idx <- which(partable$op == "~*~" & !(partable$lhs %in% lv.names) &
                     partable$block == g)
        out[idx] <- 1.0
    }

    # 5a ":="
    idx <- which(partable$op == ":=")
    if(length(idx) > 0L) {
        x <- out[ partable$free & !duplicated(partable$free) ]
        out[idx] <- lavobject@Model@def.function(x)
    }

    # 5b "=="
    idx <- which(partable$op == "==")
    if(length(idx) > 0L) {
        x <- out[ partable$free & !duplicated(partable$free) ]
        out[idx] <- lavobject@Model@ceq.function(x)
    }

    # 5c. "<" or ">"
    idx <- which((partable$op == "<" | partable$op == ">"))
    if(length(idx) > 0L) {
        x <- out[ partable$free & !duplicated(partable$free) ]
        out[idx] <- lavobject@Model@cin.function(x)
    }

    out
}


standardize.est.all.nox <- function(lavobject, partable=NULL, est=NULL,
                                    est.std=NULL, GLIST = NULL,
                                    cov.std = TRUE) {

    if(is.null(partable)) partable <- lavobject@ParTable
    if(is.null(est))   est <- lav_object_inspect_est(lavobject)
    if(is.null(est.std)) {
        est.std <- standardize.est.lv(lavobject, partable = partable, est = est,
                                      GLIST = GLIST, cov.std = cov.std)
    }
    if(is.null(GLIST)) GLIST <- lavobject@Model@GLIST

    out <- est.std; N <- length(est.std)
    stopifnot(N == length(partable$lhs))

    VY <- computeVY(lavmodel = lavobject@Model, GLIST = GLIST,
                    diagonal.only = TRUE)


    for(g in 1:lavobject@Model@nblocks) {

        ov.names     <- vnames(lavobject@ParTable, "ov",     block = g) 
        ov.names.x   <- vnames(lavobject@ParTable, "ov.x",   block = g)
        ov.names.nox <- vnames(lavobject@ParTable, "ov.nox", block = g)
        lv.names     <- vnames(lavobject@ParTable, "lv",     block = g)

        OV2 <- VY[[g]]
        # replace zero values by NA (but keep negative values)
        zero.idx <- which(abs(OV2) < .Machine$double.eps)
        if(length(zero.idx) > 0L) {
            OV2[zero.idx] <- as.numeric(NA)
        }

        # replace negative values by NA (for sqrt)
        tmp.OV2 <- OV2
        neg.idx <- which(tmp.OV2 < 0)
        if(length(neg.idx) > 0L) {
            tmp.OV2[neg.idx] <- as.numeric(NA)
        }
        OV  <- sqrt(tmp.OV2)

        if(lavobject@Model@conditional.x) {
            # extend OV with ov.names.x
            ov.names.x <- vnames(lavobject@ParTable, "ov.x", block = g)
            ov.names <- c(ov.names.nox, ov.names.x)
            OV2 <- c(OV2, diag(lavobject@SampleStats@cov.x[[g]]))
            OV <- c(OV, sqrt(diag(lavobject@SampleStats@cov.x[[g]])))
        }

        # 1a. "=~" regular indicators
        idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
                     partable$block == g)
        out[idx] <- out[idx] / OV[ match(partable$rhs[idx], ov.names) ]

        # 1b. "=~" regular higher-order lv indicators

        # 1c. "=~" indicators that are both in ov and lv
        #idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
        #                             & partable$rhs %in% lv.names &
        #             partable$block == g)

        # 2. "~" regressions (and "<~")
        idx <- which((partable$op == "~" | partable$op == "<~") & 
                     partable$lhs %in% ov.names &
                     partable$block == g)
        out[idx] <- out[idx] / OV[ match(partable$lhs[idx], ov.names) ]

        idx <- which((partable$op == "~" | partable$op == "<~") & 
                     partable$rhs %in% ov.names.nox &
                     partable$block == g)
        out[idx] <- out[idx] * OV[ match(partable$rhs[idx], ov.names.nox) ]

        # 3a. "~~" ov
        # ATTENTION: in Mplus 4.1, the off-diagonal residual covariances 
        #            were computed by the formula cov(i,j) / sqrt(i.var*j.var)
        #            were i.var and j.var where diagonal elements of OV
        #
        #            in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
        #            elements are the 'THETA' diagonal elements!!

        # variances
        rv.idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) & 
                        !(partable$lhs %in% ov.names.x) &
                        partable$lhs == partable$rhs &
                        partable$block == g)
        #out[rv.idx] <- ( out[rv.idx] / OV[ match(partable$lhs[rv.idx], ov.names) ]
         #                            / OV[ match(partable$rhs[rv.idx], ov.names) ] )
        out[rv.idx] <- ( out[rv.idx] / 
                             OV2[ match(partable$lhs[rv.idx], ov.names) ] )

        # covariances ov
        # three types:
        # - only lhs is OV (and fixed.x = FALSE)
        # - only rhs is OV (and fixed.x = FALSE)
        # - both lhs and rhs are OV (regular case)
        if(cov.std) {
            if(!is.complex(est[rv.idx])) {
                RV <- sqrt(abs(est[rv.idx]))
            } else {
                RV <- sqrt( est[rv.idx] )
            }
            rv.names <- partable$lhs[rv.idx]
        }

        # left
        idx.lhs <- which(partable$op == "~~" &
                         !(partable$lhs %in% lv.names) &
                         !(partable$lhs %in% ov.names.x) &
                         partable$lhs != partable$rhs &
                         partable$block == g)
        if(length(idx.lhs) > 0L) {
            if(cov.std == FALSE) {
                out[idx.lhs] <- 
                   (out[idx.lhs] / OV[ match(partable$lhs[idx.lhs], ov.names)])
            } else {
                out[idx.lhs] <- 
                   (out[idx.lhs] / RV[ match(partable$lhs[idx.lhs], rv.names)])
            }
        }

        # right
        idx.rhs <- which(partable$op == "~~" & 
                         !(partable$rhs %in% lv.names) &
                         !(partable$rhs %in% ov.names.x) &
                         partable$lhs != partable$rhs &
                         partable$block == g)
        if(length(idx.rhs) > 0L) {
            if(cov.std == FALSE) {
                out[idx.rhs] <- 
                    (out[idx.rhs] / OV[ match(partable$rhs[idx.rhs], ov.names)])
            } else {
                out[idx.rhs] <- 
                    (out[idx.rhs] / RV[ match(partable$rhs[idx.rhs], rv.names)])
            }
        }

        # 3b. "~~" lv
        #idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
        #             partable$block == g)

        # 4a. "~1" ov
        idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
                     !(partable$lhs %in% ov.names.x) &
                     partable$block == g)
        out[idx] <- out[idx] / OV[ match(partable$lhs[idx], ov.names) ]

        # 4b. "~1" lv
        #idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
        #             partable$block == g)

        # 4c. "|" thresholds
        idx <- which(partable$op == "|" & !(partable$lhs %in% lv.names) &
                     partable$block == g)
        out[idx] <- out[idx] / OV[ match(partable$lhs[idx], ov.names) ]

        # 4d. "~*~" scales
        idx <- which(partable$op == "~*~" & !(partable$lhs %in% lv.names) &
                     partable$block == g)
        out[idx] <- 1.0
    }

    # 5a ":="
    idx <- which(partable$op == ":=")
    if(length(idx) > 0L) {
        x <- out[ partable$free & !duplicated(partable$free) ]
        out[idx] <- lavobject@Model@def.function(x)
    }

    # 5b "=="
    idx <- which(partable$op == "==")
    if(length(idx) > 0L) {
        x <- out[ partable$free & !duplicated(partable$free) ]
        out[idx] <- lavobject@Model@ceq.function(x)
    }

    # 5c. "<" or ">"
    idx <- which((partable$op == "<" | partable$op == ">"))
    if(length(idx) > 0L) {
        x <- out[ partable$free & !duplicated(partable$free) ]
        out[idx] <- lavobject@Model@cin.function(x)
    }

    out
}

unstandardize.est.ov <- function(partable, ov.var=NULL, cov.std=TRUE) {

    # check if ustart is missing; if so, look for est
    if(is.null(partable$ustart)) 
        partable$ustart <- partable$est
  
    # check if block is missing
    if(is.null(partable$block))  {
        partable$block <- rep(1L, length(partable$ustart))
    }

    stopifnot(!any(is.na(partable$ustart)))
    est <- out <- partable$ustart
    N <- length(est)

    # nblocks
    nblocks <- lav_partable_nblocks(partable)

    # if ov.var is NOT a list, make a list
    if(!is.list(ov.var)) {
        tmp <- ov.var
        ov.var <- vector("list", length=nblocks)
        ov.var[1:nblocks] <- list(tmp)
    }

    for(g in 1:nblocks) {

        ov.names <- vnames(partable, "ov", block = g) # not user
        lv.names <- vnames(partable, "lv", block = g)

        OV  <- sqrt(ov.var[[g]])

        # 1a. "=~" regular indicators
        idx <- which(partable$op == "=~" & !(partable$rhs %in% lv.names) &
                     partable$block == g)
        out[idx] <- out[idx] * OV[ match(partable$rhs[idx], ov.names) ]

        # 1b. "=~" regular higher-order lv indicators

        # 1c. "=~" indicators that are both in ov and lv
        #idx <- which(partable$op == "=~" & partable$rhs %in% ov.names
        #                             & partable$rhs %in% lv.names &
        #             partable$block == g)

        # 2. "~" regressions (and "<~")
        idx <- which((partable$op == "~" | partable$op == "<~") & 
                     partable$lhs %in% ov.names &
                     partable$block == g)
        out[idx] <- out[idx] * OV[ match(partable$lhs[idx], ov.names) ]

        idx <- which((partable$op == "~" | partable$op == "<~") & 
                     partable$rhs %in% ov.names &
                     partable$block == g)
        out[idx] <- out[idx] / OV[ match(partable$rhs[idx], ov.names) ]

        # 3a. "~~" ov
        # ATTENTION: in Mplus 4.1, the off-diagonal residual covariances 
        #            were computed by the formula cov(i,j) / sqrt(i.var*j.var)
        #            were i.var and j.var where diagonal elements of OV
        #
        #            in Mplus 6.1 (but also AMOS and EQS), the i.var and j.var
        #            elements are the 'THETA' diagonal elements!!

        # variances
        rv.idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) & 
                        partable$lhs == partable$rhs &
                        partable$block == g)
        out[rv.idx] <- ( out[rv.idx] * OV[ match(partable$lhs[rv.idx], ov.names) ]
                                     * OV[ match(partable$rhs[rv.idx], ov.names) ] )

        # covariances
        idx <- which(partable$op == "~~" & !(partable$lhs %in% lv.names) &
                     partable$lhs != partable$rhs &
                     partable$block == g)
        if(length(idx) > 0L) {
            if(cov.std == FALSE) {
                out[idx] <- ( out[idx] * OV[ match(partable$lhs[idx], ov.names) ]
                                       * OV[ match(partable$rhs[idx], ov.names) ] )
            } else {
                if(!is.complex(out[rv.idx])) {
                    RV <- sqrt(abs(out[rv.idx]))
                } else {
                    RV <- sqrt( out[rv.idx] )
                }
                rv.names <- partable$lhs[rv.idx]
                out[idx] <- ( out[idx] * RV[ match(partable$lhs[idx], rv.names) ]
                                       * RV[ match(partable$rhs[idx], rv.names) ] )
            }
        }

        # 3b. "~~" lv
        #idx <- which(partable$op == "~~" & partable$rhs %in% lv.names &
        #             partable$block == g)

        # 4a. "~1" ov
        idx <- which(partable$op == "~1" & !(partable$lhs %in% lv.names) &
                     partable$block == g)
        out[idx] <- out[idx] * OV[ match(partable$lhs[idx], ov.names) ]

        # 4b. "~1" lv
        #idx <- which(partable$op == "~1" & partable$lhs %in% lv.names &
        #             partable$block == g)

    }

    # 5a ":="
    # 5b "=="
    # 5c. "<" or ">"

    out
}