##
##  c o m p l e x s t e p . R  Complex Step Derivation
##


complexstep <- function(f, x0, h = 1e-20, test = FALSE, ...) {
    fun <- match.fun(f)
    f <- function(x) fun(x, ...)

    f_csd <- Im(f(x0 + h * 1i)) / h
    # Test it as analycity may fail for function f
    if (test) {
        h0 <- 1e-5  # 2 * .Machine$double.eps^(1/3)
        f_grd <- (f(x0+h0/2) - f(x0-h0/2))/ h0
        if (abs(f_csd - f_grd) > h0)
            warning(paste("Complex Step check failed: ", abs(f_csd - f_grd)))
        else
            cat("Derivation Difference: ", abs(f_csd - f_grd), "\n")
    }
    return(f_csd)
}


complexstepJ <- function(f, x0, h = 1e-20, ...) {
    fun <- match.fun(f)
    f <- function(x) fun(x, ...)

    z <- f(x0)
    n <- length(x0); m <- length(z)
    J <- matrix(NA, nrow = m, ncol = n)

    for (k in 1:n) {
        x1 <- x0
        x1[k] <- x1[k] + h * 1i
        J[ , k] <- Im(f(x1)) / h
    }
    # drop matrix dimensions if n = m = 1
    if (m == 1 && n == 1) J <- J[,]
    return(J)
}