# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
"plot.qsreg" <- function(x, pch = "*", main = NA, 
    ...) {
    out <- x
    old.par <- par("mfrow", "oma")
    on.exit(par(old.par))
    set.panel(2, 2, relax = TRUE)
    plot(out$x, out$y, xlab = "X", ylab = "y", pch = pch)
    orderx <- order(out$x)
    temp <- out$fitted.values[, c(out$ind.cv, out$ind.cv.ps)]
    matlines(out$x[orderx], temp[orderx, ], lty = 1, col = c(1, 
        2))
    ##
    # residual plot
    #
    matplot(out$x, qsreg.psi(out$residuals[, c(out$ind.cv, out$ind.cv.ps)], 
        out$alpha, out$sc), col = c(1, 2), pch = "o", ylab = "Pseudo residuals", 
        xlab = "X")
    yline(0)
    if (nrow(out$cv.grid) > 1) {
        ind <- out$cv.grid[, 3] < 1e+19
        out$cv.grid <- out$cv.grid[ind, ]
        matplot(out$cv.grid[, 2], cbind(out$cv.grid[, 3], out$cv.grid[, 
            6]), xlab = "Effective number of parameters", ylab = "Log CV Rho function ", 
            log = "y", type = "l", col = c(1, 2))
        xline(out$cv.grid[out$ind.cv, 2], col = 1)
        xline(out$cv.grid[out$ind.cv.ps, 2], col = 2)
        title(" CV curves", cex = 0.5)
    }
    bplot(qsreg.psi(out$residuals[, c(out$ind.cv, out$ind.cv.ps)], 
        out$alpha, out$sc), names = c("CV", "CV pseudo"))
    yline(0, col = 2)
    if (is.na(main)) 
        mtext(deparse(out$call), cex = 1.3, outer = TRUE, line = -2)
    else mtext(main, cex = 1.3, outer = TRUE, line = -2)
}
"predict.qsreg" <- function(object, x, derivative = 0, 
    model = object$ind.cv.ps, ...) {
    if (missing(x)) 
        x <- object$x
    c(splint(object$predicted$x, object$predicted$y[, model], 
        x, derivative = derivative))
}
"qsreg" <- function(x, y, lam = NA, maxit = 50, maxit.cv = 10, 
    tol = 1e-07, offset = 0, sc = sqrt(var(y)) * 1e-05, alpha = 0.5, 
    wt = rep(1, length(x)), cost = 1, nstep.cv = 80, hmin = NA, 
    hmax = NA, trmin = 2 * 1.05, trmax = 0.95 * length(unique(x))) {
    out <- list()
    class(out) <- c("qsreg")
    N <- length(y)
    out$N <- N
    xgrid <- sort(unique(x))
    if (length(x) != length(y)) 
        stop(" X and Y do not match")
    if (!is.na(lam[1])) 
        hgrid <- log(lam)
    else {
        # find lambda grid
        if (is.na(hmin)) {
            hmin <- 0
            for (k in 1:25) {
                b <- qsreg.trace(lam = as.double(exp(hmin)), 
                  x = x, y = y, wt = wt, cost = cost, maxit = maxit, 
                  tol = tol, sc = sc, alpha = alpha)
                if (b > trmax) {
                  break
                }
                hmin <- hmin - 1
            }
        }
        if (is.na(hmax)) {
            hmax <- 0
            for (k in 1:25) {
                b <- qsreg.trace(lam = as.double(exp(hmax)), 
                  x = x, y = y, wt = wt, cost = cost, maxit = maxit, 
                  tol = tol, sc = sc, alpha = alpha)
                if (b < trmin) {
                  break
                }
                hmax <- hmax + 1
            }
        }
        h <- seq(hmin, hmax, , nstep.cv)
        lam <- exp(h)
    }
    # now loop through values for lam ( really log lam)
    b <- list()
    NL <- length(lam)
    NG <- length(xgrid)
    h <- log(lam)
    residuals <- matrix(NA, ncol = NL, nrow = N)
    diagA <- residuals
    cv.ps <- rep(0, NL)
    trace.ps <- rep(0, NL)
    cv <- rep(0, NL)
    predicted <- matrix(NA, ncol = NL, nrow = NG)
    trace <- rep(0, NL)
    converge <- rep(0, NL)
    wt.old <- wt
    for (k in 1:NL) {
        b <- .Fortran("rcss", h = as.double(h[k]), npoint = as.integer(N), 
            x = as.double(x), y = as.double(y), wt = as.double(wt.old), 
            sy = as.double(rep(0, N)), trace = as.double(0), 
            diag = as.double(rep(0, N)), cv = as.double(0), ngrid = as.integer(NG), 
            xg = as.double(xgrid), yg = as.double(rep(0, NG)), 
            job = as.integer(c(3, 3, 0)), ideriv = as.integer(0), 
            din = as.double(c(cost, offset, maxit, tol, sc, alpha)), 
            dout = as.double(rep(0, 4)), ierr = as.integer(0))
        residuals[, k] <- y - b$sy
        diagA[, k] <- b$diag
        cv[k] <- b$dout[4]
        trace[k] <- b$trace
        predicted[, k] <- b$yg
        converge[k] <- b$dout[1]
        wt.old <- b$wt
    }
    # second loop to find approx CV residuals based on pseudo values
    y.pseudo <- rep(NA, N)
    residuals.cv <- matrix(NA, ncol = NL, nrow = length(x))
    for (k in 1:NL) {
        y.pseudo <- (sc) * qsreg.psi(residuals[, k], alpha = alpha, 
            C = sc) + y - residuals[, k]
        #
        # call the robust spline but set the cutoff for the huber weight so big
        # it is essentially a LS spline this helps to match the lambda for robust spline
        # with a LS one.
        #
        b <- .Fortran("rcss", h = as.double(h[k]), npoint = as.integer(N), 
            x = as.double(x), y = as.double(y.pseudo), wt = as.double(wt), 
            sy = as.double(rep(0, N)), trace = as.double(0), 
            diag = as.double(rep(0, N)), cv = as.double(0), ngrid = as.integer(NG), 
            xg = as.double(xgrid), yg = as.double(rep(0, NG)), 
            job = as.integer(c(3, 3, 0)), ideriv = as.integer(0), 
            din = as.double(c(cost, offset, maxit, tol, sqrt(var(y)) * 
                10, alpha)), dout = as.double(rep(0, 4)), ierr = as.integer(0))

        #
        # CV residuals based on pseudo-data)
        # Use the linear approximation  Y_k - f.cv_k = (Y_k- f_k)/( 1- A_kk)
        #  f.cv_k = f_k/( 1- A_kk) - ( A_kk)Y_k/( 1- A_kk)
        #
        # Note: we find f.cv based on pseudo data but then consider its deviation
        # from the actual data
        #
        f.cv <- (b$sy/(1 - b$diag)) - b$diag * y.pseudo/(1 - 
            b$diag)
        trace.ps[k] <- b$trace
        residuals.cv[, k] <- (y - f.cv)
        cv.ps[k] <- mean(qsreg.rho(y - f.cv, alpha = alpha, C = sc))
    }
    #
    #
    #
    cv.grid <- cbind(lam, trace, cv, converge, trace.ps, cv.ps)
    dimnames(cv.grid) <- list(NULL, c("lambda", "trace", "CV", 
        "iterations", "trace.PS", "CV.PS"))
    #
    ind.cv <- (1:NL)[cv == min(cv)]
    ind.cv.ps <- (1:NL)[cv.ps == min(cv.ps)]
    out$call <- match.call()
    out$x <- x
    out$y <- y
    out$predicted <- list(x = xgrid, y = predicted)
    out$trace <- trace
    out$residuals.cv <- residuals.cv
    out$residuals <- residuals
    out$fitted.values <- y - residuals
    out$cv.grid <- cv.grid
    out$diagA <- diagA
    out$sc <- sc
    out$alpha <- alpha
    out$ind.cv <- ind.cv
    out$ind.cv.ps <- ind.cv.ps
    out
}
"qsreg.fit" <- function(x, y, lam, maxit = 50, maxit.cv = 10, 
    tol = 1e-04, offset = 0, sc = sqrt(var(y)) * 1e-07, alpha = 0.5, 
    wt = rep(1, length(x)), cost = 1) {
    N <- length(y)
    if (length(x) != length(y)) 
        stop(" X and Y do not match")
    h <- log(lam)
    temp <- .Fortran("rcss", h = as.double(log(lam)), npoint = as.integer(N), 
        x = as.double(x), y = as.double(y), wt = as.double(wt), 
        sy = as.double(rep(0, N)), trace = as.double(0), diag = as.double(rep(0, 
            N)), cv = as.double(0), ngrid = as.integer(0), xg = as.double(0), 
        yg = as.double(0), job = as.integer(c(3, 0, 0)), ideriv = as.integer(0), 
        din = as.double(c(cost, offset, maxit, tol, sc, alpha)), 
        dout = as.double(rep(0, 4)), ierr = as.integer(0))$dout
    return(temp)
}
"qsreg.psi" <- function(r, alpha = 0.5, C = 1) {
    temp <- rep(NA, length(r))
    r <- r/C
    temp <- r
    ind <- r > 1
    temp[ind] <- 2 * alpha
    ind <- r < 1 & r > 0
    temp[ind] <- (2 * alpha * r[ind])
    ind <- r < -1
    temp[ind] <- -2 * (1 - alpha)
    ind <- r > -1 & r < 0
    temp[ind] <- 2 * (1 - alpha) * r[ind]
    temp
}
"qsreg.rho" <- function(r, alpha = 0.5, C = 1) {
    temp <- rep(NA, length(r))
    ind <- r > C
    temp[ind] <- 2 * alpha * r[ind] - alpha * C
    ind <- r < C & r > 0
    temp[ind] <- (alpha * r[ind]^2)/C
    ind <- r < -C
    temp[ind] <- -2 * (1 - alpha) * r[ind] - (1 - alpha) * C
    ind <- r > -C & r < 0
    temp[ind] <- ((1 - alpha) * r[ind]^2)/C
    temp
}