R/reparametrize.R

Defines functions .reparametrize.cov .reparametrize.atanh .reparametrize.logvar .reparametrize.var .reparametrize.logsd .reparametrize.sd .reparametrize.logsquare .reparametrize.square .reparametrize.log .reparametrize.one .reparametrize.none .init_transform .reparametrize reparametrize

### reparametrize.R --- 
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: Apr 25 2021 (11:22) 
## Version: 
## Last-Updated: jun 16 2022 (15:03) 
##           By: Brice Ozenne
##     Update #: 735
##----------------------------------------------------------------------
## 
### Commentary: 
## 
### Change Log:
##----------------------------------------------------------------------
## 
### Code:

## * reparametrize
reparametrize <- function(p, type, level, sigma, k.x, k.y, 
                          Jacobian = TRUE, dJacobian = TRUE,
                          FUN = NULL,
                          transform.sigma = NULL,
                          transform.k = NULL,
                          transform.rho = NULL,
                          transform.names = TRUE){

    n.p <- length(p)
    name.p <- names(p)
    options <- LMMstar.options()
        
    ## ** normalize user arguments
    if(!is.null(FUN)){
        if(is.function(FUN)==FALSE){
            stop("Argument \'FUN\' must either be null or a function. \n")
        }
        if(!is.null(transform.sigma)){
            warning("Argument \'transform.sigma\' is ignored when argument \'FUN\' is specified. \n")
        }
        if(!is.null(transform.k)){
            warning("Argument \'transform.k\' is ignored when argument \'FUN\' is specified. \n")
        }
        if(!is.null(transform.rho)){
            warning("Argument \'transform.rho\' is ignored when argument \'FUN\' is specified. \n")
        }
        valid.args <- c("inverse","p","type","sigma","k.x","k.y")

        if(any(names(formals(FUN)) %in% c(valid.args,"...") == FALSE)){
            stop("Invalid argument for \'FUN\'. \n",
                 "Invalid argument(s): \'",paste(names(formals(FUN))[names(formals(FUN)) %in% valid.args == FALSE], collapse = "\' \'"),"\' \n",
                 "Valid arguments: \'",paste(valid.args, collapse = "\' \'"),"\' \n")
        }
        if(any(valid.args %in% names(formals(FUN)) == FALSE) && "..." %in% names(formals(FUN)) == FALSE){
            stop("Missing argument for \'FUN\': ",paste(valid.args[valid.args %in% names(formals(FUN)) == FALSE], collapse = "\' \'"),"\' \n")
        }
        if(as.numeric(dJacobian) %in% 0:2 == FALSE){
            stop("Argument \'dJacobian\' must be 0 (FALSE), 1 (TRUE), or 2. \n")
        }
        
        ## *** transformed parameter
        p.trans <- FUN(p = p, type = type, sigma = sigma, k.x = k.x, k.y = k.y, inverse = FALSE)
        newname <- names(p.trans)
        names(p.trans) <- names(p)
        if(length(p.trans)!=n.p){
            stop("Argument \'FUN\' must output a vector of length ",n.p," (same length as argument \'p\'). \n")
        }

        if(transform.names){
            attr(p.trans,"newname") <- newname
        }

        ## *** jacobian
        if(Jacobian){
            if(any(abs(p-FUN(p = p.trans, type = type, sigma = sigma, k.x = k.x, k.y = k.y, inverse = TRUE))>1e-6)){
                stop("The argument \'inverse\' of the function defined in argument \'FUN\' seems to lead to incorrect value. \n",
                     "FUN(FUN(p, inverse = FALSE),inverse = TRUE) should be equal to argument \'p\'. \n")
            }

            dFUN <- function(x, type, sigma, k.x, k.y, vectorize){
                out <- numDeriv::jacobian(function(pp){FUN(p = pp, type = type, sigma = sigma, k.x = k.x, k.y = k.y, inverse = TRUE)},x)
                if(vectorize){
                    return(as.vector(out))
                }else{
                    return(out)
                }
            }
            attr(p.trans,"Jacobian") <- dFUN(p.trans, type = type, sigma = sigma, k.x = k.x, k.y = k.y, vectorize = FALSE)
            dimnames(attr(p.trans,"Jacobian")) <- list(name.p,name.p)
        }
        
        ## *** derivative of the jacobian
        if(dJacobian==1){
            ddFUN <- function(x, type, sigma, k.x, k.y, vectorize){
                out <- numDeriv::jacobian(function(pp){dFUN(x = pp, type = type, sigma = sigma, k.x = k.x, k.y = k.y, vectorize = TRUE)},x)
                if(vectorize){ ## convert to array
                    return(out)
                }else{
                    A.out <- array(NA, dim = rep(n.p, 3), dimnames = list(name.p,name.p,name.p))
                    for(iCol in 1:NCOL(out)){
                        A.out[,,iCol] <- out[,iCol]
                    }
                    return(A.out)
                }
            }        
            attr(p.trans,"dJacobian") <- ddFUN(p.trans, type = type, sigma = sigma, k.x = k.x, k.y = k.y, vectorize = FALSE)
        }else if(dJacobian==2){
            stop("Argument \'dJacobian==2\' not implemented for numerical derivatives. \n")
        }

    }else{
        init <- .init_transform(transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho, 
                                x.transform.sigma = options$transform.sigma, x.transform.k = options$transform.k, x.transform.rho = options$transform.rho)
        attr(init$transform.sigma,"arg") <- transform.sigma
        attr(init$transform.k,"arg") <- transform.k
        attr(init$transform.rho,"arg") <- transform.rho

        ls.out <- .reparametrize(p = p, type = type, level = level, sigma = sigma, k.x = k.x, k.y = k.y,
                                 Jacobian = Jacobian, dJacobian = dJacobian, inverse = FALSE,
                                 transform.sigma = init$transform.sigma,
                                 transform.k = init$transform.k,
                                 transform.rho = init$transform.rho,
                                 transform.names = transform.names)

        p.trans <- ls.out$p
        if(transform.names){
            attr(p.trans,"newname") <- ls.out$newname
        }
        if(Jacobian){
            attr(p.trans,"Jacobian") <- ls.out$Jacobian
        }
        if(dJacobian>0){
            attr(p.trans,"dJacobian") <- ls.out$dJacobian
        }
    }
    
    ## ** export
    return(p.trans)
}

## * .reparametrize
.reparametrize <- function(p, type, level,
                           sigma, k.x, k.y,
                           Jacobian, dJacobian, inverse,
                           transform.sigma,
                           transform.k,
                           transform.rho,
                           transform.names){

    arg.transform.sigma <- attr(transform.sigma,"arg")
    arg.transform.k <- attr(transform.k,"arg")

    if(transform.rho %in% "cov" && any("rho" %in% type)){
        transform.sigma <- "square"
        transform.k <- "var"
        if(!is.null(arg.transform.sigma) && arg.transform.sigma!=transform.sigma){
            warning("Argument \'transform.sigma\' ignored when argument \'transform.rho\' is \"cov\". \n")
        }
        if(!is.null(arg.transform.k) && arg.transform.k!=transform.k){
            warning("Argument \'transform.k\' ignored when argument \'transform.rho\' is \"cov\". \n")
        }
    }
    if(transform.k %in% c("sd","logsd","var","logvar") && any("k" %in% type)){
        transform.sigma <- switch(transform.k,
                                  "sd" = "none",
                                  "logsd" = "log",
                                  "var" = "square",
                                  "logvar" = "logsquare")

        if(!is.null(arg.transform.sigma) && arg.transform.sigma!=transform.sigma){
            warning(paste0("Argument \'transform.sigma\' ignored when argument \'transform.k\' is \"",transform.k,"\". \n"))
        }
    }

    ## ** initialize
    n.p <- length(p)
    name.p <- names(p)
    index.sigma <- which(type=="sigma")
    index.k <- which(type=="k")
    index.rho <- which(type=="rho")
    out <- list(p = p,
                newname = name.p,
                transform = (transform.sigma!="none") || (transform.k!="none") || (transform.rho!="none"),
                transform.sigma = transform.sigma,
                transform.k = transform.k,
                transform.rho = transform.rho)
    if(Jacobian){
        out$Jacobian <- matrix(0, nrow = n.p, ncol = n.p, dimnames = list(name.p,name.p))
        diag(out$Jacobian) <- 1
    }        
    if(dJacobian>0){
        out$dJacobian <- array(0, dim = rep(n.p, 3), dimnames = list(name.p,name.p,name.p))
    }
    
 
    ## *** sigma (or k with log/square/logsquare)
    if(length(index.sigma)>0){

        transform.sigma <- switch(transform.sigma,
                                  "none" = .reparametrize.none,
                                  "one" = .reparametrize.one,
                                  "log" = .reparametrize.log,
                                  "square" = .reparametrize.square,
                                  "logsquare" = .reparametrize.logsquare)

        if(!transform.names){
            level.sigma <- NULL
        }else if(transform.k %in% c("sd","logsd","var","logvar")){
            level.sigma <- sapply(level,"[[",1)
        }else{
            level.sigma <- stats::setNames(rep(NA,n.p),name.p)
            pos.column <- which(grepl(":",name.p,fixed=TRUE))
            if(length(pos.column)>1){ ## add strata
                level.sigma[pos.column] <- sapply(strsplit(name.p[pos.column],":", fixed = TRUE), function(iString){
                    paste(c(":",iString[-1]),collapse="")
                })
            }
        }
        out <- transform.sigma(p = p, out = out,
                               index = index.sigma,
                               level = level.sigma, type = "sigma",
                               inverse = inverse, transform.names = transform.names, Jacobian = Jacobian, dJacobian = dJacobian)
               
    }

    ## *** k
    if(length(index.k)>0){

        transform.k <- switch(transform.k,
                              "none" = .reparametrize.none,
                              "log" = .reparametrize.log,
                              "square" = .reparametrize.square,
                              "logsquare" = .reparametrize.logsquare,
                              "sd" = .reparametrize.sd,
                              "logsd" = .reparametrize.logsd,
                              "var" = .reparametrize.var,
                              "logvar" = .reparametrize.logvar)

        if(!transform.names){
            level.k <- NULL
        }else{
            level.k <- sapply(level,"[[",1)
        }
        out <- transform.k(p = p, out = out,
                           index = index.k, indexSigma = match(sigma[index.k], name.p),
                           level = level.k, type = "k",
                           inverse = inverse, transform.names = transform.names, Jacobian = Jacobian, dJacobian = dJacobian)
               

    }

    ## *** rho
    if(length(index.rho)>0){

        transform.rho <- switch(transform.rho,
                                "none" = .reparametrize.none,
                                "atanh" = .reparametrize.atanh,
                                "cov" = .reparametrize.cov)

        if(!transform.names){
            level.rho <- NULL
        }else{
            level.rho <- sapply(level,"[[",1)
        }
        out <- transform.rho(p = p, out = out,
                             index = index.rho, indexSigma = match(sigma[index.rho], name.p), indexKx = match(k.x[index.rho], name.p), indexKy = match(k.y[index.rho], name.p),
                             level = level.rho, type = "rho",
                             inverse = inverse, transform.names = transform.names, Jacobian = Jacobian, dJacobian = dJacobian)

    }     

    ## ** export
    return(out)
}

## * .init_transform
.init_transform <- function(transform.sigma, transform.k, transform.rho, 
                            x.transform.sigma, x.transform.k, x.transform.rho){

    ## ** transform
    if(is.null(transform.sigma)){
        transform.sigma.save <- transform.sigma
        transform.sigma <- x.transform.sigma
        attr(transform.sigma,"arg") <- transform.sigma.save
    }else{
        if(identical(transform.sigma,"")){
            transform.sigma <- "none"
        }
        transform.sigma <- match.arg(transform.sigma, c("none","one","log","square","logsquare"))
        attr(transform.sigma,"arg") <- NULL
    }
    if(is.null(transform.k)){
        transform.k.save <- transform.k
        transform.k <- x.transform.k
        attr(transform.k,"arg") <- transform.k.save
    }else{
        if(identical(transform.k,"")){
            transform.k <- "none"
        }
        transform.k <- match.arg(transform.k, c("none","log","square","logsquare","sd","logsd","var","logvar"))
        attr(transform.k,"arg") <- NULL
    }
    if(is.null(transform.rho)){
        transform.rho.save <- transform.rho
        transform.rho <- x.transform.rho
        attr(transform.rho,"arg") <- transform.rho.save
    }else{
        if(identical(transform.rho,"")){
            transform.rho <- "none"
        }
        transform.rho <- match.arg(transform.rho, c("none","atanh", "cov"))
        attr(transform.rho,"arg") <- NULL
    }

    ## ** x.transform
    if(!is.null(x.transform.sigma) && !is.null(x.transform.k) && !is.null(x.transform.rho)){
        if(is.function(transform.sigma) || is.function(transform.k) || is.function(transform.rho)){
            test.notransform <- FALSE
        }else{
            test.notransform <- (transform.sigma==x.transform.sigma) && (transform.k==x.transform.k) && (transform.rho==x.transform.rho)
        }
    }else{
        test.notransform <- NULL
    }

    ## ** export
    return(list(transform.sigma = transform.sigma,
                transform.k = transform.k,
                transform.rho = transform.rho,
                test.notransform = test.notransform))
}

## * .reparametrize.none
.reparametrize.none <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                                transform.names, Jacobian, dJacobian){

    if(!inverse && transform.names && all(!is.na(level[index]))){
        out$newname[index] <- paste0(type, level[index])
    }

    return(out)

}

## * .reparametrize.one
.reparametrize.one <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                               transform.names, Jacobian, dJacobian){

    out$p[index] <- 1

    if(transform.names){
        out$newname[index] <- "reference"
    }

    if(Jacobian){
        diag(out$Jacobian)[index] <- 0
    }

    return(out)

}

## * .reparametrize.log
.reparametrize.log <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                               transform.names, Jacobian, dJacobian){

    if(inverse){

        ## f^{-1}:y -> x=exp(y)
        out$p[index] <- exp(p[index])

    }else{
           
        ## f:x -> y=log(x)
        out$p[index] <- log(p[index])
            
        if(transform.names){
            if(all(is.na(level[index]))){
                out$newname[index] <- paste0("log(",type,")")
            }else{
                out$newname[index] <- paste0("log(",type,")", level[index])
            }
        }
            
        if(Jacobian){
            ## df^{-1}:y -> exp(y)=x
            diag(out$Jacobian)[index] <- p[index] 
        }

        if(dJacobian==1){ ## orignal scale after tranformation
            for(iIndex in index){
                ## d2f^{-1}:y -> exp(y)=x
                out$dJacobian[iIndex,iIndex,iIndex] <- p[iIndex]
            }
        }else if(dJacobian==2){ ## newscale after transformation
            for(iIndex in index){
                ## d/dx df^{-1}:y -> 1
                out$dJacobian[iIndex,iIndex,iIndex] <- 1
            }
        }
    }

    return(out)

}

## * .reparametrize.square
.reparametrize.square <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                                  transform.names, Jacobian, dJacobian){
    if(inverse){

        ## f^{-1}:y -> x=sqrt(y)
        out$p[index] <- sqrt(p[index])

    }else{

        ## f:x -> y=x^2
        out$p[index] <- p[index]^2
                
        if(transform.names){
            if(all(is.na(level[index]))){
                out$newname[index] <- paste0(type, "^2")
            }else{
                out$newname[index] <- paste0(type, "^2", level[index])
            }
        }

        if(Jacobian){
            ## df^{-1}:y -> 1/(2*sqrt(y))=1/(2x)
            diag(out$Jacobian)[index] <- 1/(2*p[index]) 
        }

        if(dJacobian == 1){
            for(iIndex in index){
                ## d2f^{-1}:y -> -1/(4*y^(3/2))=-1/(4x^3)
                out$dJacobian[iIndex,iIndex,iIndex] <- -1/(4*p[iIndex]^3)
            }
        }else if(dJacobian == 2){
            for(iIndex in index){
                ## d/dx df^{-1}:y -> -1/(2x^2)
                out$dJacobian[iIndex,iIndex,iIndex] <- -1/(2*p[iIndex]^2)
            }
        }
    }

    return(out)
}

## * .reparametrize.logsquare
.reparametrize.logsquare <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                                     transform.names, Jacobian, dJacobian){
    if(inverse){

        ## f^{-1}:y -> x=sqrt(exp(y))=exp(y/2)
        out$p[index] <- exp(p[index]/2)

    }else{

        ## f:x -> y=log(x^2)=2log(x)
        out$p[index] <- log(p[index]^2)

        if(transform.names){
            if(all(is.na(level[index]))){
                out$newname[index] <- paste0("log(",type,"^2)")
            }else{
                out$newname[index] <- paste0("log(",type,"^2)", level[index])
            }
        }

        if(Jacobian){
            ## df^{-1}:y -> x=d exp(y/2)= exp(y/2)/2 = x/2
            diag(out$Jacobian)[index] <- p[index]/2 
        }

        if(dJacobian == 1){
            for(iIndex in index){
                ## d2f^{-1}:y -> exp(y/2)/4 = x/4
                out$dJacobian[iIndex,iIndex,iIndex] <- p[iIndex]/4
            }
        }else if(dJacobian == 2){
            for(iIndex in index){
                ## d/dx df^{-1}:y -> 1/2
                out$dJacobian[iIndex,iIndex,iIndex] <- 1/2
            }
        }
    }

    return(out)
}

## * .reparametrize.sd
.reparametrize.sd <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                              transform.names, Jacobian, dJacobian){

    if(type!="k"){
        stop("Only implemented for type \'k\' parameters. \n")
    }
    
    n.index <- length(index)

    if(inverse){
        ## f^{-1}:(y1,y2) -> (x1,x2)=(y1,y2/y1)
        out$p[index] <- p[index]/p[indexSigma]

    }else{

        ## f:(x1,x2) -> (y1,y2)=(x1,x1*x2)
        out$p[index] <- p[index]*p[indexSigma]
         
        if(transform.names){
            out$newname[index] <- paste0("sigma", level[index])
        }

        if(Jacobian){
            ## df^{-1}:(y1,y2) -> (1,-y2/y1^2)=  (1,-x2/x1)            /y1 
            ##                 -> (0,1/y1)    =  (0,1/x1)              /y2
            ## Here we are only interested in the second column, i.e. (-x2/x1,1/x1)
            for(iK in 1:n.index){
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$Jacobian[iIndex.k,iIndex.sigma] <- -p[iIndex.k]/p[iIndex.sigma]
            }
            diag(out$Jacobian)[index] <- 1/p[indexSigma]
        }
    }
    
    if(dJacobian == 1){
        for(iK in 1:n.index){
            ## df^{-1}:(y1,y2) -> (-y2/y1^2,1/y1)     
            ## d2f^{-1}:(y1,y2) -> (2 y2/y1^3,-1/y1^2) = (2x2/x1^2,-1/x1^1)    
            ##                  -> (-1/y1^2,0)         = (-1/x^2,0)
            iIndex.k <- index[iK]
            iIndex.sigma <- indexSigma[iK]
            out$dJacobian[iIndex.k,iIndex.sigma,iIndex.sigma] <- 2*p[iIndex.k]/(p[iIndex.sigma]^2)
            out$dJacobian[iIndex.k,iIndex.k,iIndex.sigma] <- -1/p[iIndex.sigma]^2
            out$dJacobian[iIndex.k,iIndex.sigma,iIndex.k] <- -1/p[iIndex.sigma]^2
        }
    }else if(dJacobian == 2){
        for(iK in 1:n.index){
            ## df^{-1}:(y1,y2) -> (-x2/x1,1/x1)     
            ## d/dx df^{-1}:(y1,y2) -> (x2/x1^2,-1/x1^2)
            ##                      -> (-1/x1,0)
            iIndex.k <- index[iK]
            iIndex.sigma <- indexSigma[iK]
            out$dJacobian[iIndex.k,iIndex.sigma,iIndex.sigma] <- p[iIndex.k]/p[iIndex.sigma]^2
            out$dJacobian[iIndex.k,iIndex.k,iIndex.sigma] <- -1/p[iIndex.sigma]^2
            out$dJacobian[iIndex.k,iIndex.sigma,iIndex.k] <- -1/p[iIndex.sigma]
        }
    }

    return(out)
}

## * .reparametrize.logsd
.reparametrize.logsd <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                                 transform.names, Jacobian, dJacobian){

    if(type!="k"){
        stop("Only implemented for type \'k\' parameters. \n")
    }
    
    n.index <- length(index)

    if(inverse){

        ## f^{-1}:(y1,y2) -> (x1,x2)=(exp(y1),exp(y2-y1))
        out$p[index] <- exp(p[index]-p[indexSigma])

    }else{

        ## f:(x1,x2) -> (y1,y2)=(log(x1),log(x1*x2))
        out$p[index] <- log(p[index]*p[indexSigma])

        if(transform.names){
            out$newname[index] <- paste0("log(sigma)", level[index])
        }

        if(Jacobian){
            ## df^{-1}:(y1,y2) -> (exp(y1),-exp(y2-y1)) = (x1, -x2)            /y1
            ##                 -> (0,exp(y2-y1)) = (0, x2)                     /y2
            ## Here we are only interested in the second column, i.e. (-x2,x2)
            for(iK in 1:n.index){
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$Jacobian[iIndex.k,iIndex.sigma] <- -p[iIndex.k]
            }
            diag(out$Jacobian)[index] <- p[index]
        }
        if(dJacobian == 1){
            for(iK in 1:n.index){
                ## df^{-1}:(y1,y2) -> (-exp(y2-y1),exp(y2-y1))
                ## d2f^{-1}:(y1,y2) -> (exp(y2-y1),-exp(y2-y1)) = (x2,-x2)
                ##                  -> (-exp(y2-y1),exp(y2-y1)) = (-x2,x2)
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.sigma] <- p[iIndex.k]
                out$dJacobian[iIndex.k,iIndex.k,iIndex.sigma] <- -p[iIndex.k]
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.k] <- -p[iIndex.k]
                out$dJacobian[iIndex.k,iIndex.k,iIndex.k] <- p[iIndex.k]
            }
        }else if(dJacobian == 2){
            for(iK in 1:n.index){
                ## df^{-1}:(y1,y2) -> (-x2,x2)
                ## d/dx df^{-1}:(y1,y2) -> (0,0)
                ##                      -> (-1,1)
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.k] <- -1
                out$dJacobian[iIndex.k,iIndex.k,iIndex.k] <- 1
            }                            
        }

    }

    return(out)
}

## * .reparametrize.var
.reparametrize.var <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                               transform.names, Jacobian, dJacobian){

    if(type!="k"){
        stop("Only implemented for type \'k\' parameters. \n")
    }
    
    n.index <- length(index)

    if(inverse){

        ## f^{-1}:(y1,y2) -> (x1,x2)=(sqrt(y1),sqrt(y2/y1))
        out$p[index] <- sqrt(p[index]/p[indexSigma])

    }else{

        ## f:(x1,x2) -> (y1,y2)=(x1^2,x1^2*x2^2)
        out$p[index] <- p[indexSigma]^2*p[index]^2

        if(transform.names){
            out$newname[index] <- paste0("sigma^2", level[index])
        }

        if(Jacobian){
            ## df^{-1}:(y1,y2) -> (1/(2*sqrt(y1)),-sqrt(y2)/(2*y1^(3/2))=  (1/(2*x1),-x2/(2*x1^2))            /y1 
            ##                 -> (0,1/(2*sqrt(y2*y1)))              =  (0,1/(2*x2*x1^2))                     /y2
            ## Here we are only interested in the second column, i.e. (-x2/(2*x1^2),1/(2*x1^2*x2))
            for(iK in 1:n.index){
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$Jacobian[iIndex.k,iIndex.sigma] <- -p[iIndex.k]/(2*p[iIndex.sigma]^2)
            }
            diag(out$Jacobian)[index] <- 1/(2*p[indexSigma]^2*p[index])
            
        }
        if(dJacobian == 1){
            ## sqrt(\sigma1), sqrt(\sigma2/\sigma1), ...   so  -1/(4 \sigma1^(3/2)) 3\sqrt(\sigma2)/(4 \sigma1^(5/2))      0 -1/(4 \sqrt(\sigma2) \sigma1^(3/2))        
            ##                                                      0              - 1 / (4*\sqrt(\sigma2)\sigma1^(3/2))   0 -1 / (4*\sigma2^(3/2)\sqrt(\sigma1))
            ## WARNING we plug-in sqrt(\sigma1) and sqrt(\sigma2/\sigma1)
            for(iK in 1:n.index){
                ## df^{-1}:(y1,y2) -> (-sqrt(y2)/(2*y1^(3/2)),1/(2*sqrt(y2*y1)))
                ## d2f^{-1}:(y1,y2) -> (3sqrt(y2)/(4*y1^(5/2)),-1/(4*sqrt(y2)*y1^(3/2)))) = (3 x2/(4x1^4), -1/(4*x2*x1^4))
                ##                  -> (-1/(4*sqrt(y2)*y1^(3/2)),-1/(4*y2^(3/2)*y1))) = (-1/(4*x2*x1^4), -1/(4*x2^3*x1^4))
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.sigma] <- 3*p[iIndex.k]/(4*p[iIndex.sigma]^4)
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.k] <- -1/(4*p[iIndex.sigma]^4*p[iIndex.k])
                out$dJacobian[iIndex.k,iIndex.k,iIndex.sigma] <- -1/(4*p[iIndex.sigma]^4*p[iIndex.k])
                out$dJacobian[iIndex.k,iIndex.k,iIndex.k] <- -1/(4*p[iIndex.sigma]^4*p[iIndex.k]^3)
            }
        }else if(dJacobian == 2){
            for(iK in 1:n.index){
                ## df^{-1}:(y1,y2) -> (-x2/(2*x1^2),1/(2*x1^2*x2))
                ## d/dx df^{-1}:(y1,y2) -> (x2/(4*x1^3),-1/(4*x1^3*x2))
                ##                      -> (-1/(2*x1^2),-1/(2*x1^2*x2^2))
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.sigma] <- p[iIndex.k]/(4*p[iIndex.sigma]^3)
                out$dJacobian[iIndex.k,iIndex.k,iIndex.sigma] <- -1/(4*p[iIndex.sigma]^3*p[iIndex.k])
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.k] <- -1/(2*p[iIndex.sigma]^2)
                out$dJacobian[iIndex.k,iIndex.k,iIndex.k] <- -1/(2*p[iIndex.sigma]^2*p[iIndex.k]^2)
            }
        }
    }

    return(out)
}

## * .reparametrize.logvar
.reparametrize.logvar <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                                  transform.names, Jacobian, dJacobian){
    if(type!="k"){
        stop("Only implemented for type \'k\' parameters. \n")
    }
    
    n.index <- length(index)

    if(inverse){

        ## f^{-1}:(y1,y2) -> (x1,x2)=(sqrt(exp(y1)),sqrt(exp(y2-y1)))=(exp(y1/2),exp(y2/2-y1/2))
        out$p[index] <- exp(p[index]/2-p[indexSigma]/2)

    }else{

        ## f:(x1,x2) -> (y1,y2)=(log(x1^2),log(x1^2*x2^2))=(2*log(x1),2*log(x1*x2))
        out$p[index] <- 2*log(p[index]*p[indexSigma])

        if(transform.names){
            out$newname[index] <- paste0("log(sigma^2)", level[index])
        }
        
        if(Jacobian){
            ## df^{-1}:(y1,y2) -> (exp(y1/2)/2,-exp(y2/2-y1/2)/2) = (x1/2, -x2/2)            /y1
            ##                 -> (0,exp(y2/2-y1/2)/2) = (0, x2*2)                          /y2
            ## Here we are only interested in the second column, i.e. (-x2/2/,x2/2)
            for(iK in 1:n.index){
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$Jacobian[iIndex.k,iIndex.sigma] <- -p[iIndex.k]/2
            }
            diag(out$Jacobian)[index] <- p[index]/2
        }
        if(dJacobian == 1){
            ## df^{-1}:(y1,y2) -> (-exp(y2/2-y1/2)/2,exp(y2/2-y1/2)/2)
            ## d2f^{-1}:(y1,y2) -> (exp(y2/2-y1/2)/4,-exp(y2/2-y1/2)/4) = (x2/4,-x2/4)
            ##                  -> (-exp(y2/2-y1/2)/4,exp(y2/2-y1/2)/4) = (-x2/4,x2/4)
            for(iK in 1:n.index){
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.sigma] <- p[iIndex.k]/4
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.k] <- -p[iIndex.k]/4
                out$dJacobian[iIndex.k,iIndex.k,iIndex.sigma] <- -p[iIndex.k]/4
                out$dJacobian[iIndex.k,iIndex.k,iIndex.k] <- p[iIndex.k]/4
            }
        }else if(dJacobian == 2){
            for(iK in 1:n.index){
                ## df^{-1}:(y1,y2) -> (-x2/2/,x2/2)
                ## d/dx df^{-1}:(y1,y2) -> (0,0)
                ##                      -> (-1/2,1/2)
                iIndex.k <- index[iK]
                iIndex.sigma <- indexSigma[iK]
                out$dJacobian[iIndex.k,iIndex.sigma,iIndex.k] <- -1/2
                out$dJacobian[iIndex.k,iIndex.k,iIndex.k] <- 1/2
            }            
        }

    }
    return(out)                
}

## * .reparametrize.atanh
.reparametrize.atanh <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                                 transform.names, Jacobian, dJacobian){
    if(type!="rho"){
        stop("Only implemented for type \'rho\' parameters. \n")
    }
    
    n.index <- length(index)

    if(inverse){

        ## f^{-1}:y -> x=tanh(y)
        out$p[index] <- tanh(p[index])

    }else{

        ## f:x -> y=atanh(x)
        out$p[index] <- atanh(p[index])

        if(transform.names){
            out$newname[index] <- paste0("atanh(rho)", level[index])
        }

        if(Jacobian){
            ## df^{-1}:y -> d tanh(y) = 1-tanh(y)^2 = 1-x^2
            diag(out$Jacobian)[index] <- (1-p[index]^2) ## (1-tanh(x)^2) where x=atanh(rho)
        }
        if(dJacobian == 1){
            for(iRho in index){
                ## d2f^{-1}:y -> d tanh(y) = - 2 tanh(y)(1-tanh(y)^2) = -2 x (1-x^2)
                out$dJacobian[iRho,iRho,iRho] <- -2*p[iRho]*(1-p[iRho]^2)
            }
        }else if(dJacobian == 2){
            for(iRho in index){
                ## d/dx df^{-1}:y -> -2x
                out$dJacobian[iRho,iRho,iRho] <- -2*p[iRho]
            }
        }
    }    
    
    return(out)
}

## * .reparametrize.cov
.reparametrize.cov <- function(p, out, index, indexSigma, indexKx, indexKy, level, inverse, type,
                               transform.names, Jacobian, dJacobian){

    if(type!="rho"){
        stop("Only implemented for type \'rho\' parameters. \n")
    }
    
    n.index <- length(index)

    Sigma <- p[indexSigma]
    Kx <- p[indexKx]
    Ky <- p[indexKy]
        
    if(inverse){

        ## fa^{-1}:(y1,y4) -> (x1,x4)=(sqrt(y1),y4/y1)
        ## fb^{-1}:(y1,y2,y3) -> (x1,x2,x3)=(sqrt(y1),sqrt(y2/y1),y4/sqrt(y1*y2))
        ## fc^{-1}:(y1,y2,y3) -> (x1,x2,x3)=(sqrt(y1),sqrt(y2/y1),y4/y2)
        ## fd^{-1}:(y1,y2,y3,y4) -> (x1,x2,x3,x4)=(sqrt(y1),sqrt(y2/y1),sqrt(y3/y1),y4/sqrt(y2*y3))
        Kx[is.na(Kx)] <- Sigma[is.na(Kx)]
        Ky[is.na(Ky)] <- Sigma[is.na(Ky)]
        out$p[index] <- p[index]/sqrt(Kx*Ky)

    }else{

        ## fa:(x1,x4) -> (y1,y4)=(x1^2,x1^2*x4)
        ## fb:(x1,x2,x4) -> (y1,y2,y4)=(x1^2,x1^2*x2^2,x1^2*x2*x4)
        ## fc:(x1,x2,x4) -> (y1,y2,y4)=(x1^2,x1^2*x2^2,x1^2*x2^2*x4)
        ## fd:(x1,x2,x3,x4) -> (y1,y2,y3,y4)=(x1^2,x1^2*x2^2,x1^2*x3^2,x1^2*x2*x3*x4)
        Kx[is.na(Kx)] <- 1
        Ky[is.na(Ky)] <- 1
        out$p[index] <- p[index]*Sigma^2*Kx*Ky

        if(transform.names){
            out$newname[index] <- paste0("cov", level[index])
        }

        if(Jacobian){
            ## dfa^{-1}:(y1,y4) -> (1/(2*sqrt(y1)),-y4/y1^2)=  (1/(2*x1),-x4/x1^2)          /y1 
            ##                  -> (0,1/y1)                 =  (0,1/x1^2)                   /y4
            ## Here we are only interested in the second column, i.e. (-x4/(x1^2),1/(x1^2))

            ## dfb^{-1}:(y1,y2,y4) -> -y4/(2*sqrt(y2)*y1^(3/2)) =  -x4 /(2*y1) = -x4 / (2*x1^2)      /y1 
            ##                     -> -y4/(2*sqrt(y1)*y2^(3/2)) =  -x4 /(2*y2) = -x4 / (2*x1^2*x2^2) /y2
            ##                     -> 1/sqrt(y1*y2)             =  1/(x1^2*x2)                       /y4

            ## dfc^{-1}:(y1,y2,y4) -> 0                                        /y1 
            ##                     -> -y4/y2^2 = -x4/y2 = -x4/(x1^2*x2^2)      /y2
            ##                     -> 1/y2              =  1/(x1^2*x2^2)       /y4

            ## dfd^{-1}:(y1,y2,y3,y4) -> 0                                                            /y1 
            ##                        -> -y4/(2*sqrt(y3)*y2^(3/2)) = -x4/(2*y2) = -x4/(2*x1^2*x2^2)   /y2
            ##                        -> -y4/(2*sqrt(y2)*y3^(3/2)) = -x4/(2*y3) = -x4/(2*x1^2*x3^2)   /y3
            ##                        -> 1/sqrt(y2*y3)     =  1/(x1^2*x2*x3)                          /y4

            for(iRho in 1:n.index){ ## iRho <- 1
                iIndex.sigma <- indexSigma[iRho]
                iIndex.kx <- indexKx[iRho]
                iIndex.ky <- indexKy[iRho]
                iIndex.rho <- index[iRho]
                if(is.na(iIndex.kx) && is.na(iIndex.ky)){
                    out$Jacobian[iIndex.rho,iIndex.sigma] <- -p[iIndex.rho]/p[iIndex.sigma]^2 ## /y1
                }else if(!is.na(iIndex.kx) && is.na(iIndex.ky)){
                    out$Jacobian[iIndex.rho,iIndex.sigma] <- -p[iIndex.rho]/(2*p[iIndex.sigma]^2) ## /y1                    
                    out$Jacobian[iIndex.rho,iIndex.kx] <- -p[iIndex.rho]/(2*p[iIndex.sigma]^2*p[iIndex.kx]^2) ## /y2                    
                }else if(is.na(iIndex.kx) && !is.na(iIndex.ky)){
                    out$Jacobian[iIndex.rho,iIndex.sigma] <- -p[iIndex.rho]/(2*p[iIndex.sigma]^2) ## /y1                    
                    out$Jacobian[iIndex.rho,iIndex.ky] <- -p[iIndex.rho]/(2*p[iIndex.sigma]^2*p[iIndex.ky]^2) ## /y2                   
                }else{ ## both non-NA
                    ## 0 for /y1
                    out$Jacobian[iIndex.rho,iIndex.kx] <- out$Jacobian[iIndex.rho,iIndex.sigma]-p[iIndex.rho]/(2*p[iIndex.sigma]^2*p[iIndex.kx]^2) ## /y3
                    out$Jacobian[iIndex.rho,iIndex.ky] <- out$Jacobian[iIndex.rho,iIndex.sigma]-p[iIndex.rho]/(2*p[iIndex.sigma]^2*p[iIndex.ky]^2) ## /y4                   
                }
            }
            diag(out$Jacobian)[index] <- 1/(Sigma^2*Kx*Ky) ## /y4 ok
            
        }
        if(dJacobian == 1){
            ##    fa^{-1}:(y1,y4) -> (x1,x4)=(sqrt(y1),y4/y1)
            ##    dfa^{-1}:(y1,y4) -> (-y4/y1^2,1/y1)
            ## so d2fa^{-1}:(y1,y4) -> (2 y4/y1^3,-1/y1^2)   =  (2 x4/x1^4,-1/x1^4)     // y1
            ##                      -> (-1/y1^2,0)           = (-1/x1^4,0)             // y4

            ##    fb^{-1}:(y1,y2,y3) -> (x1,x2,x3)=(sqrt(y1),sqrt(y2/y1),y4/sqrt(y1*y2))
            ##    dfb^{-1}:(y1,y2,y4) -> ( -y4/(2*sqrt(y2)*y1^(3/2)),-y4/(2*sqrt(y1)*y2^(3/2)),1/sqrt(y1*y2) )
            ## so d2fb^{-1}:(y1,y2,y4) -> (3*y4/(4*sqrt(y2)*y1^(5/2)), y4/(4*y1^(3/2)*y2^(3/2)), -1/(2*y1^(3/2)*sqrt(y2))) = (3*x4/(4*x1^4)   ,x4/(4*x1^4*x2^2) ,-1/(2*x1^4*x2)) ## y1
            ##                         -> (y4/(4*y1^(3/2)*y2^(3/2)), 3*y4/(4*sqrt(y1)*y2^(5/2)), -1/(2*sqrt(y1)*y2^(3/2))) = (x4/(4*x1^4*x2^2),3*x4/(4*x1^4*x2^4),-1/(2*x1^4*x2^3)) ## y2
            ##                         -> (-1/(2*sqrt(y2)*y1^(3/2)),-1/(2*sqrt(y1)*y2^(3/2)) , 0)                          = (-1/(2*x1^4*x2)  , -1/(2*x1^4*x2^3), 0)     ## y4
        
            ## fd^{-1}:(y1,y2,y3,y4) -> (x1,x2,x3,x4)=(sqrt(y1),sqrt(y2/y1),sqrt(y3/y1),y4/sqrt(y2*y3))
            ## dfd^{-1}:(y1,y2,y3,y4) -> (0, -y4/(2*sqrt(y3)*y2^(3/2)), -y4/(2*sqrt(y2)*y3^(3/2)), 1/sqrt(y2*y3))
            ## so d2fd^{-1}:(y1,y2,y3,y4) -> (0, 0, 0, 0)
            ##                            -> (0, 3*y4/(4*sqrt(y3)*y2^(5/2)), y4/(4*y2^(3/2)*y3^(3/2)), -1/(2y2^(3/2)*sqrt(y3))) = (0,3*x4/(4*x2^4*x1^4),x4/(4*x2^2*x3^2*x1^4),-1/(2*x2^3*x3*x1^4)))
            ##                            -> (0, y4/(4*y3^(3/2)*y2^(3/2)), 3y4/(4*sqrt(y2)*y3^(5/2)), -1/(2sqrt(y2)*y3^(3/2))) = (0,x4/(4*x3^2^x2^2*x1^4),3*x4/(4*y3^4*x1^4),-1/(2*x2*x3^3*x1^4))
            ##                            -> (0, -1/(2*sqrt(y3)*y2^(3/2)), -1/(2*sqrt(y2)*y3^(3/2)), 0) = (0,-1/(2*x3*x2^3*x1^4),-1/(2*x2*x3^3*x1^4),0)

            for(iRho in 1:n.index){ ## iRho <- 1
                iIndex.sigma <- indexSigma[iRho]
                iIndex.kx <- indexKx[iRho]
                iIndex.ky <- indexKy[iRho]
                iIndex.rho <- index[iRho]
                if(is.na(iIndex.kx) && is.na(iIndex.ky)){

                    out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.sigma] <- 2*p[iIndex.rho]/p[iIndex.sigma]^4
                    out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.rho] <- -1/p[iIndex.sigma]^4
                    out$dJacobian[iIndex.rho,iIndex.rho,iIndex.sigma] <- -1/p[iIndex.sigma]^4

                }else if(!is.na(iIndex.kx) && is.na(iIndex.ky)){
                    
                    out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.sigma] <- 3*p[iIndex.rho]/(4*p[iIndex.sigma]^4)               
                    out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.kx] <- out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.kx] + p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.kx]^2)               
                    out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.rho] <- -1/(2*p[iIndex.sigma]^4*p[iIndex.kx])               

                    out$dJacobian[iIndex.rho,iIndex.kx,iIndex.sigma] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.sigma] + p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.kx]^2)
                    out$dJacobian[iIndex.rho,iIndex.kx,iIndex.kx] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.kx] + 3*p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.kx]^4)
                    out$dJacobian[iIndex.rho,iIndex.kx,iIndex.rho] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.rho] - 1/(2*p[iIndex.sigma]^4*p[iIndex.kx]^3)

                    out$dJacobian[iIndex.rho,iIndex.rho,iIndex.sigma] <- -1/(2*p[iIndex.sigma]^4*p[iIndex.kx])
                    out$dJacobian[iIndex.rho,iIndex.rho,iIndex.kx] <- out$dJacobian[iIndex.rho,iIndex.rho,iIndex.kx] - 1/(2*p[iIndex.sigma]^4*p[iIndex.kx]^3)


                }else if(is.na(iIndex.kx) && !is.na(iIndex.ky)){

                    out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.sigma] <- 3*p[iIndex.rho]/(4*p[iIndex.sigma]^4)               
                    out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.ky] + p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.ky]^2)               
                    out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.rho] <- -1/(2*p[iIndex.sigma]^4*p[iIndex.ky])               

                    out$dJacobian[iIndex.rho,iIndex.ky,iIndex.sigma] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.sigma] + p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.ky]^2)
                    out$dJacobian[iIndex.rho,iIndex.ky,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.ky] + 3*p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.ky]^4)
                    out$dJacobian[iIndex.rho,iIndex.ky,iIndex.rho] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.rho] - 1/(2*p[iIndex.sigma]^4*p[iIndex.ky]^3)

                    out$dJacobian[iIndex.rho,iIndex.rho,iIndex.sigma] <- -1/(2*p[iIndex.sigma]^4*p[iIndex.ky])
                    out$dJacobian[iIndex.rho,iIndex.rho,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.rho,iIndex.ky] - 1/(2*p[iIndex.sigma]^4*p[iIndex.ky]^3)

                }else{ ## both non-NA
                    ## 0 for /y1
                    out$dJacobian[iIndex.rho,iIndex.kx,iIndex.kx] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.kx]+3*p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.kx]^4)
                    out$dJacobian[iIndex.rho,iIndex.kx,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.ky]+p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.kx]^2*p[iIndex.ky]^2)
                    out$dJacobian[iIndex.rho,iIndex.kx,iIndex.rho] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.rho]-1/(2*p[iIndex.sigma]^4*p[iIndex.kx]^3*p[iIndex.ky])

                    out$dJacobian[iIndex.rho,iIndex.ky,iIndex.kx] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.kx]+p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.kx]^2*p[iIndex.ky]^2)
                    out$dJacobian[iIndex.rho,iIndex.ky,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.ky]+3*p[iIndex.rho]/(4*p[iIndex.sigma]^4*p[iIndex.ky]^4)
                    out$dJacobian[iIndex.rho,iIndex.ky,iIndex.rho] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.rho]-1/(2*p[iIndex.sigma]^4*p[iIndex.kx]*p[iIndex.ky]^3)

                    out$dJacobian[iIndex.rho,iIndex.rho,iIndex.kx] <- out$dJacobian[iIndex.rho,iIndex.rho,iIndex.kx]-1/(2*p[iIndex.sigma]^4*p[iIndex.kx]^3*p[iIndex.ky]) 
                    out$dJacobian[iIndex.rho,iIndex.rho,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.rho,iIndex.ky]-1/(2*p[iIndex.sigma]^4*p[iIndex.kx]*p[iIndex.ky]^3) 
                }
                
            }
            
            
        }else if(dJacobian == 2){

            for(iRho in 1:n.index){
                ## dfa^{-1}:(y1,y4) -> (-x4/(x1^2),1/(x1^2))
                ## d/dx dfa^{-1}:(y1,y4) -> (2*x4/(x1^3),-2/(x1^3))
                ##                       -> (-1/(x1^2),0)

                ## dfb^{-1}:(y1,y2,y4) -> (-x4 / (2*x1^2), -x4 / (2*x1^2*x2^2), 1/(x1^2*x2))
                ## d/dx dfb^{-1}:(y1,y2,y4) -> (x4 / (x1^3),   x4 / (x1^3*x2^2),   -2/(x1^3*x2))
                ##                          -> (0,             x4 / (x1^2*x2^3),   -1/(x1^2*x2^2))
                ##                          -> (-1 / (2*x1^2), -1 / (2*x1^2*x2^2),  0)
                
                ## dfd^{-1}:(y1,y2,y3,y4) -> (0, -x4/(2*x1^2*x2^2), -x4/(2*x1^2*x3^2), 1/(x1^2*x2*x3))
                ## d/dx dfd^{-1}:(y1,y2,y3,y4) -> (0, x4/(x1^3*x2^2),   x4/(x1^3*x3^2),  -2/(x1^3*x2*x3))
                ##                             -> (0, x4/(x1^2*x2^3),   0,               -1/(x1^2*x2^2*x3))
                ##                             -> (0, 0,                x4/(x1^2*x3^3),  -1/(x1^2*x2*x3^2))
                ##                             -> (0, -1/(2*x1^2*x2^2), -1/(2*x1^2*x3^2), 0)
                for(iRho in 1:n.index){ ## iRho <- 1
                    iIndex.sigma <- indexSigma[iRho]
                    iIndex.kx <- indexKx[iRho]
                    iIndex.ky <- indexKy[iRho]
                    iIndex.rho <- index[iRho]
                    if(is.na(iIndex.kx) && is.na(iIndex.ky)){
                        out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.sigma] <- 2*p[iIndex.rho]/(p[iIndex.sigma]^3)
                        out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.rho] <- -1/p[iIndex.sigma]^2
                    }else if(!is.na(iIndex.kx) && is.na(iIndex.ky)){
                        out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.sigma] <- p[iIndex.rho]/(p[iIndex.sigma]^3)
                        out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.rho] <- -1/(2*p[iIndex.sigma]^2)
                    
                        out$dJacobian[iIndex.rho,iIndex.kx,iIndex.sigma] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.sigma]+p[iIndex.rho]/(p[iIndex.sigma]^3*p[iIndex.kx]^2)
                        out$dJacobian[iIndex.rho,iIndex.kx,iIndex.kx] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.kx]+p[iIndex.rho]/(p[iIndex.sigma]^2*p[iIndex.kx]^3)
                        out$dJacobian[iIndex.rho,iIndex.kx,iIndex.rho] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.rho]-1/(2*p[iIndex.sigma]^2*p[iIndex.kx]^2)
                    }else if(is.na(iIndex.kx) && !is.na(iIndex.ky)){
                        out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.sigma] <- p[iIndex.rho]/(p[iIndex.sigma]^3)
                        out$dJacobian[iIndex.rho,iIndex.sigma,iIndex.rho] <- -1/(p[iIndex.sigma]^2)
                   
                        out$dJacobian[iIndex.rho,iIndex.ky,iIndex.sigma] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.sigma]+p[iIndex.rho]/(p[iIndex.sigma]^3*p[iIndex.ky]^2)
                        out$dJacobian[iIndex.rho,iIndex.ky,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.ky]+p[iIndex.rho]/(p[iIndex.sigma]^2*p[iIndex.ky]^3)
                        out$dJacobian[iIndex.rho,iIndex.ky,iIndex.rho] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.rho]-1/(2*p[iIndex.sigma]^2*p[iIndex.ky]^2)
                    }else{ ## both non-NA
                        out$dJacobian[iIndex.rho,iIndex.kx,iIndex.sigma] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.sigma]+p[iIndex.rho]/(p[iIndex.sigma]^3*p[iIndex.kx]^2) 
                        out$dJacobian[iIndex.rho,iIndex.kx,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.ky]+p[iIndex.rho]/(p[iIndex.sigma]^2*p[iIndex.kx]^3) 
                        out$dJacobian[iIndex.rho,iIndex.kx,iIndex.rho] <- out$dJacobian[iIndex.rho,iIndex.kx,iIndex.rho]-1/(2*p[iIndex.sigma]^2*p[iIndex.kx]^2) 
                        
                        out$dJacobian[iIndex.rho,iIndex.ky,iIndex.sigma] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.sigma]+p[iIndex.rho]/(p[iIndex.sigma]^3*p[iIndex.ky]^2)
                        out$dJacobian[iIndex.rho,iIndex.ky,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.ky]+p[iIndex.rho]/(p[iIndex.sigma]^2*p[iIndex.ky]^3)
                        out$dJacobian[iIndex.rho,iIndex.ky,iIndex.rho] <- out$dJacobian[iIndex.rho,iIndex.ky,iIndex.rho]-1/(2*p[iIndex.sigma]^2*p[iIndex.ky]^2)
                    }
                    
                    out$dJacobian[iIndex.rho,iIndex.rho,iIndex.sigma] <- -2/(p[iIndex.sigma]^3*p[iIndex.kx]*p[iIndex.ky])
                    if(!is.na(iIndex.kx)){
                        out$dJacobian[iIndex.rho,iIndex.rho,iIndex.kx] <- out$dJacobian[iIndex.rho,iIndex.rho,iIndex.kx]-1/(p[iIndex.sigma]^2*p[iIndex.kx]^2*p[iIndex.ky])
                    }
                    if(!is.na(iIndex.ky)){
                        out$dJacobian[iIndex.rho,iIndex.rho,iIndex.ky] <- out$dJacobian[iIndex.rho,iIndex.rho,iIndex.ky]-1/(p[iIndex.sigma]^2*p[iIndex.kx]*p[iIndex.ky]^2)
                    }
                    

                    
                }
            }
        }
    }
    
    return(out)
}

##----------------------------------------------------------------------
### reparametrize.R ends here

Try the LMMstar package in your browser

Any scripts or data that you put into this service are public.

LMMstar documentation built on Nov. 9, 2023, 1:06 a.m.