psindex: Parameter Stability Index

# numerical derivatives using complex numbers
# see Squire & Trapp 1998, siam rev 40(1) 110-112
# or Ridout, MS (2009), the american statistician 63(1) 66-74

# it would seem that you can choose h to be fairly small, without
# sacrifycing accuracy due to rounding errors

# YR 17 July 2012

lav_func_gradient_complex <- function(func, x, 
                                      h = .Machine$double.eps, ... , 
                                      check.scalar = TRUE, 
                                      fallback.simple = TRUE) {

    # check current point, see if it is a scalar function
    if(check.scalar) {
        f0 <- try(func(x*(0+1i), ...), silent = TRUE)
        if(inherits(f0, "try-error")) {
            if(fallback.simple) {
                dx <- lav_func_gradient_simple(func = func, x = x, h = sqrt(h),
                                               check.scalar = check.scalar, ...)
                return(dx)
            } else {
                stop("function does not support non-numeric (complex) argument")
            }
        }
        if(length(f0) != 1L) {
            stop("function is not scalar and returns more than one element")
        }
    }

    nvar <- length(x)

    # determine 'h' per element of x
    h <- pmax(h, abs(h*x))

    # get exact h, per x
    tmp <- x + h
    h <- (tmp - x)

    # simple 'forward' method
    dx <- rep(as.numeric(NA), nvar)
    for(p in seq_len(nvar)) {
        dx[p] <- Im(func(x + h*1i*(seq.int(nvar) == p),...))/h[p]
    }

    dx
}

# as a backup, if func() is not happy about non-numeric arguments
lav_func_gradient_simple <- function(func, x, 
                                     h = sqrt(.Machine$double.eps), ... , 
                                     check.scalar = TRUE) {

    # check current point, see if it is a scalar function
    if(check.scalar) {
        f0 <- func(x, ...)
        if(length(f0) != 1L) {
            stop("function is not scalar and returns more than one element")
        }
    }

    nvar <- length(x)

    # determine 'h' per element of x
    h <- pmax(h, abs(h*x))

    # get exact h, per x
    tmp <- x + h
    h <- (tmp - x)

    # simple 'forward' method
    dx <- rep(as.numeric(NA), nvar)
    for(p in seq_len(nvar)) {
        dx[p] <- (func(x + h*(seq.int(nvar) == p), ...) - func(x,...))/h[p]
    }

    dx
}

lav_func_jacobian_complex <- function(func, x, 
                                      h = .Machine$double.eps, ... ,
                                      fallback.simple = TRUE) {

    f0 <- try(func(x*(0+1i), ...), silent = TRUE)
    if(inherits(f0, "try-error")) {
        if(fallback.simple) {
            dx <- lav_func_jacobian_simple(func = func, x = x, h = sqrt(h), ...)
            return(dx)
        } else {
            stop("function does not support non-numeric (complex) argument")
        }
    }
    nres <- length(f0)
    nvar <- length(x)

    # determine 'h' per element of x
    h <- pmax(h, abs(h*x))

    # get exact h, per x
    tmp <- x + h
    h <- (tmp - x)

    # simple 'forward' method
    dx <- matrix(as.numeric(NA), nres, nvar)
    for(p in seq_len(nvar)) {
        dx[,p] <- Im(func(x + h*1i*(seq.int(nvar) == p), ...))/h[p]
    }

    dx
}

lav_func_jacobian_simple <- function(func, x, 
                                     h = sqrt(.Machine$double.eps), ...) {

    f0 <- func(x, ...)
    nres <- length(f0)
    nvar <- length(x)

    # determine 'h' per element of x
    h <- pmax(h, abs(h*x))

    # get exact h, per x
    tmp <- x + h
    h <- (tmp - x)

    # simple 'forward' method
    dx <- matrix(as.numeric(NA), nres, nvar)
    for(p in seq_len(nvar)) {
        dx[,p] <- (func(x + h*(seq.int(nvar) == p), ...) - func(x,...))/h[p]
    }

    dx
}

# this is based on the Ridout (2009) paper, and the code snippet for 'h4'
lav_func_hessian_complex <- function(func, x, 
                                     h = .Machine$double.eps, ... , 
                                     check.scalar = TRUE) {

    # check current point, see if it is a scalar function
    if(check.scalar) {
        f0 <- try(func(x*(0+1i), ...), silent = TRUE)
        if(inherits(f0, "try-error")) {
            stop("function does not support non-numeric (complex) argument")
        }
        if(length(f0) != 1L) {
            stop("function is not scalar and returns more than one element")
        }
    }

    nvar <- length(x)

    # determine 'h' per element of x
    #delta1 <- pmax(h^(1/3), abs(h^(1/3)*x))
    #delta2 <- pmax(h^(1/5), abs(h^(1/5)*x))
    delta1 <- h^(1/3)
    delta2 <- h^(1/5)

    H <- matrix(as.numeric(NA), nvar, nvar)
    for(i in seq_len(nvar)) {
        for(j in 1:i) {
            if(i == j) {
                delta <- delta2
            } else {
                delta <- delta1
            }
            H[i,j] <- H[j,i] <-
                Im(func(x + delta*1i*(seq.int(nvar) == i)*x + 
                            delta*(seq.int(nvar) == j)*x, ...) -
                   func(x + delta*1i*(seq.int(nvar) == i)*x - 
                            delta*(seq.int(nvar) == j)*x, ...)) / 
                   (2*delta*delta*x[i]*x[j])
        }
    }

    H
}

# quick and dirty (FIXME!!!) way to get
# surely there must be a more elegant way?
# dCor/dCov
lav_deriv_cov2cor <- function(COV = NULL, num.idx = NULL) {

    # dCor/dvar1 = - cov / (2*var1 * sqrt(var1) * sqrt(var2))
    # dCor/dvar2 = - cov / (2*var2 * sqrt(var1) * sqrt(var2))
    # dCor/dcov  =  1/(sqrt(var1) * sqrt(var2))

    # diagonal: diag(lav_matrix_vech(tcrossprod(1/delta)))

    nvar <- ncol(COV);  pstar <- nvar*(nvar+1)/2
    delta <- sqrt(diag(COV))
    if(length(num.idx) > 0L) {
        delta[num.idx] <- 1.0
    }

    A <- COV * -1/( 2*delta*delta*tcrossprod(delta) )
    if(length(num.idx) > 0L) {
        A[num.idx,] <- 0; A[cbind(num.idx, num.idx)] <- 1
    }
    A2 <- diag(nvar) %x% t(A)

    OUT <- diag( pstar )
    diag(OUT) <- lav_matrix_vech(tcrossprod(1/delta))
    var.idx <- lav_matrix_diagh_idx(nvar)
    DUP <- lav_matrix_duplication(nvar)
    OUT[,var.idx] <- t(DUP) %*% A2[,lav_matrix_diag_idx(nvar)]

    if(length(num.idx) > 0L) {
        var.idx <- var.idx[-num.idx]
    }
    OUT[var.idx, var.idx] <- 0

    OUT
}


lav_deriv_cov2cor_numerical <- function(COV, num.idx=integer(0)) {

    compute.R <- function(x) {
        S <- lav_matrix_vech_reverse(x)
        diagS <- diag(S); delta <- 1/sqrt(diagS)
        if(length(num.idx) > 0L) {
            delta[num.idx] <- 1.0
        }
        R <- diag(delta) %*% S %*% diag(delta)
        #R <- cov2cor(S)
        R.vec <- lav_matrix_vech(R, diagonal = TRUE)
        R.vec
    }

    x <- lav_matrix_vech(COV, diagonal = TRUE)
    dx <- lav_func_jacobian_complex(func=compute.R, x=x)

    dx
}