# fields, Tools for spatial data
# Copyright 2004-2007, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
"Krig.engine.default" <- function(out, verbose = FALSE) {
    #
    # matrix decompositions for computing estimate
    #
    # Computational outline:( '.' is used for subscript)
    #
    # The form of the estimate is
    #    fhat(x) = sum phi.j(x) d.j  + sum psi.k(x) c.k
    #
    # the {phi.j} are the fixed part of the model usually low order polynomials
    # and is also referred to as spatial drift.
    #
    # the {psi.k} are the covariance functions evaluated at the unique observation
    # locations or 'knots'.  If xM.k is the kth unique location psi.k(x)= k(x, xM.k)
    # xM is also out$knots in the code below.
    #
    # the goal is find decompositions that facilitate rapid solution for
    # the vectors d and c. The eigen approach below was identified by
    # Wahba, Bates Wendelberger and is stable even for near colinear covariance
    # matrices.
    # This function does the main computations leading to the matrix decompositions.
    # With these decompositions the coefficients of the solution are found in
    # Krig.coef and the GCV and REML functions in Krig.gcv.
    #
    #  First is an outline calculations with equal weights
    #  T the fixed effects regression matrix  T.ij = phi.j(xM.i)
    #  K the covariance matrix for the unique locations
    # From the spline literature the solution solves the well known system
    # of two eqautions:
    #    -K( yM - Td - Kc) + lambda *Kc = 0
    #                 -T^t ( yM-Td -Kc) = 0
    #
    # Mulitple through by K inverse and substitute, these are equivalent to
    #
    #  -1-   -( yM- Td - Kc) + lambda c = 0
    #  -2-                        T^t c = 0
    #
    #
    #  A QR decomposition is done for   T= (Q.1,Q.2)R
    #   by definition  Q.2^T T =0
    #
    #  equation  -2- can be thought of as a constraint
    # with  c= Q.2 beta2
    # substitute in  -1-  and multiply through by Q.2^T
    #
    #      -Q.2^T yM  + Q.2^T K Q.2 beta2  + lambda beta2 = 0
    #
    #   Solving
    #   beta2 = {Q.2^T K Q.2 + lambda I )^ {-1} Q.2^T yM
    #
    # and so one sloves this linear system for beta2 and then uses
    #     c= Q.2 beta2
    #   to determine c.
    #
    #  eigenvalues and eigenvectors are found for M= Q.2^T K Q.2
    #     M = V diag(eta) V^T
    #  and these facilitate solving this system efficiently for
    #  many different values of lambda.
    #  create eigenvectors, D = (0, 1/eta)
    #  and G= ( 0,0) %*% diag(D)
    #         ( 0,V)
    # so that
    #
    #          beta2 = G%*% ( 1/( 1+ lambda D)) %*% u
    # with
    #
    #          u = (0, V Q.2^T W2 yM)
    #
    # Throughout keep in mind that M has smaller dimension than G due to
    # handling the null space.
    #
    # Now solve for d.
    #
    # From -1-  Td = yM - Kc - lambda c
    #      (Q.1^T) Td =  (Q.1^T) ( yM- Kc)
    #
    #   ( lambda c is zero by -2-)
    #
    #   so Rd = (Q.1^T) ( yM- Kc)
    # use qr functions to solve triangular system in R to find d.
    #
    #----------------------------------------------------------------------
    # What about errors with a general precision matrix, W?
    #
    # This is an important case because with replicated observations the
    # problem will simplify into a smoothing problem with the replicate group
    # means and unequal measurement error variances.
    #
    # the equations to solve are
    #     -KW( yM - Td - Kc) + lambda *Kc = 0
    #     -T^t W( yM-Td -Kc) =0
    #
    # Multiple through by K inverse and substitute, these are equivalent to
    #
    #  -1b-      -W( yM- Td - Kc) + lambda c = 0
    #  -2b-      (WT)^t c = 0
    #
    # Let W2 be the symmetric square root of W,  W= W2%*% W2
    # and W2.i be the inverse of W2.
    #
    #  -1c-      -(  W2 yM - W2 T d - (W2 K W2) W2.ic) + lambda W2.i c = 0
    #  -2c-      (W2T)^t  W2c = 0
    Tmatrix <- do.call(out$null.function.name, c(out$null.args, 
        list(x = out$xM, Z = out$ZM)))
    if (verbose) {
        cat(" Model Matrix: spatial drift and Z", fill = TRUE)
        print(Tmatrix)
    }
    # Tmatrix premultiplied by sqrt of wieghts
    Tmatrix <- out$W2 %d*% Tmatrix
    qr.T <- qr(Tmatrix)
    #
    #verbose block
    if (verbose) {
        cat("first 5 rows of qr.T$qr", fill = TRUE)
        print(qr.T$qr[1:5, ])
    }
    #
    # find  Q_2 K Q_2^T  where K is the covariance matrix at the knot points
    #
    tempM <- t(out$W2 %d*% do.call(out$cov.function.name, c(out$args, 
        list(x1 = out$knots, x2 = out$knots))))
    tempM <- out$W2 %d*% tempM
    tempM <- qr.yq2(qr.T, tempM)
    tempM <- qr.q2ty(qr.T, tempM)
    np <- nrow(out$knots)
    nt <- (qr.T$rank)
    if (verbose) {
        cat("np, nt", np, nt, fill = TRUE)
    }
    #
    # Full set of decompositions for
    # estimator for nonzero lambda
    tempM <- eigen(tempM, symmetric = TRUE)
    D <- c(rep(0, nt), 1/tempM$values)
    #
    # verbose block
    if (verbose) {
        cat("eigen values:", fill = TRUE)
        print(D)
    }
    #
    # Find the transformed data vector used to
    # evaluate the solution, GCV, REML  at different lambdas
    #
    u <- c(rep(0, nt), t(tempM$vectors) %*% qr.q2ty(qr.T, c(out$W2 %d*% 
        out$yM)))
    #
    #
    return(list(D = D, qr.T = qr.T, decomp = "WBW", V = tempM$vectors, 
        u = u, nt = nt, np = np))
}