# 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

mKrig<- function (x,y,
                  weights = rep( 1, nrow( x)), lambda=0,
                  cov.function= "stationary.cov", m=2,
                  chol.args = NULL, cov.args=NULL,...)
{

# grab other arguments for covariance function
    cov.args<- c( cov.args, list( ...))
   

# see comments in Krig.engine.fixed for algorithmic commentary
#    
#
# default values for Cholesky decomposition, these are important
# for sparse matrix decompositions 

    if( is.null( chol.args)) {
        chol.args<- list( pivot= FALSE)}
    else{
         chol.args<- chol.args}
    
# check for duplicate x's.
# stop if there are any
        if( any(duplicated( cat.matrix(x))) ) 
                  stop("locations are not unique see help(mKrig) ")

# create fixed part of model as m-1 order polynomial 
         Tmatrix <- fields.mkpoly(x, m)
        
# set some dimensions
         np <- nrow( x)
         nt <- ncol( Tmatrix)

# covariance matrix at observation locations
# NOTE: if cov.function is a sparse constuct then tempM will be sparse.
# see e.g. wendland.cov

         tempM<-  do.call(
             cov.function, c(cov.args, list(x1 = x, x2 = x))  )

# record sparsity of tempM 

        if( !is.matrix( tempM)) {  
             nzero<- length( tempM@entries) }
         else{
              nzero<- np^2}

# add diagonal matrix that is the observation error Variance 
# NOTE: diag should be an overloaded function for  tempM sparse.

      if( lambda != 0) {
        diag(tempM) <-   (lambda/weights) + diag(tempM)
    
      }

# cholesky decoposition of tempM
# do.call used to supply other arguments to the function
# especially for sparse applications.
# If chol.args is NULL then this is the same as
#              Mc<-chol(tempM)
    
     Mc<-  do.call("chol",c( list( x=tempM), chol.args) ) 

# Efficent way to multply inverse of Mc times the Tmatrix
    
 VT<- forwardsolve( Mc, x=Tmatrix, transpose=TRUE, upper.tri=TRUE)
 qr( VT)-> qr.VT

# start linear algebra to find solution 
# Note that all these expressions make sense if y is a matrix
# of several data sets and one is solving for the coefficients 
# of all of these at once. In this case d.coef and c.coef are matrices

#
# now do generalized least squares for d 

 d.coef<- qr.coef( qr.VT,
                  forwardsolve( Mc, transpose=TRUE, y, upper.tri=TRUE) )

# and then find c. 
# find the coefficents for the spatial part.    
 
 c.coef<- forwardsolve( Mc, transpose=TRUE,
                        y - Tmatrix%*% d.coef, upper.tri=TRUE)
 c.coef<- backsolve( Mc,c.coef)
    

# return coefficients and   include lambda as a check because 
# results are meaningless for other values of lambda 
# returned list is an "object" of class mKrig (micro Krig)
    
    out<-  list( d=(d.coef), c=(c.coef), 
                nt=nt, 
                np=np,
                lambda.fixed=lambda,
                x= x,
                cov.function.name=cov.function, args=cov.args,m=m,
                chol.args=chol.args,
                call= match.call(), 
                nonzero.entries= nzero)
#       
   out$residuals<-  y - predict.mKrig( out)

class(out)<- "mKrig"    

return(out)

  }

print.mKrig<- function(x,...){

    c1 <- "Number of Observations:"
    c2 <- length(x$residuals)
  
        c1 <- c(c1, "Degree of polynomial null space ( base model):")
        c2 <- c(c2, x$m - 1)

    c1 <- c(c1, "Number of parameters in the null space")
    c2 <- c(c2, x$nt)

    c1 <- c(c1, "Smoothing parameter")
    c2 <- c(c2, x$lambda.fixed)

    c1 <- c(c1, "Nonzero entries in covariance")
    c2 <- c(c2, x$nonzero.entries)
     
    sum <- cbind(c1, c2)
    dimnames(sum) <- list(rep("", dim(sum)[1]), rep("", dim(sum)[2]))
    cat("Call:\n")
    dput(x$call)
    print(sum, quote = FALSE)
    cat("Covariance Model:", x$cov.function, fill = TRUE)


    if (x$cov.function == "stationary.cov") {
        cat("  Covariance function is ", x$args$Covariance, fill = TRUE)
    }
    if (!is.null(x$args)) {
        cat("  Names of non-default covariance arguments: ", 
            fill = TRUE)
        cat("      ", paste(as.character(names(x$args)), collapse = ", "), 
            fill = TRUE)
    }

    invisible(x)
  }    




predict.mKrig<- function(object, xnew=NULL, derivative=0,...){


# the main reason to pass new args to the covariance is to increase
# the temp space size for sparse multiplications
# other optional arguments from mKrig are passed along in the
# list object$args

   cov.args<- list(...)

# predict at observation locations by default
  
  if( is.null( xnew)){
        xnew<- object$x }

# fixed part of the model this a polynomial of degree m-1    
# Tmatrix <- fields.mkpoly(xnew, m=object$m)
# 

  if( derivative==0){
  temp1 <- fields.mkpoly(xnew, m=object$m) %*% object$d}
  else{
       temp1<- fields.derivative.poly( xnew,m=object$m, object$d)}

# add nonparametric part. Covariance basis functions
# times coefficients.
# syntax is the name of the function and then a list with
# all the arguments. This allows for different covariance functions
# that have been passed as their name.


  if( derivative == 0){
# argument list are the parameters and other options from mKrig
#  locations and coefficients,
  temp2 <-  do.call(object$cov.function.name, 
             c(object$args,
             list(x1 = xnew, x2 = object$x, C = object$c), cov.args ) )}
  else{
  temp2 <-  do.call(object$cov.function.name, 
             c(
                object$args,
                list(x1 = xnew, x2 = object$x, C = object$c,
                           derivative=derivative), 
                cov.args )
            )}



# add two parts together and coerce to vector 
 return( (temp1 + temp2) )

}