MLE.objective.fn <- function(ltheta, info, value = TRUE) {
    # marginal process covariance matrix
    y <- as.matrix(info$y)
    x <- info$x
    smoothness <- info$smoothness
    ngrid <- info$ngrid
    # number of reps
    M <- ncol(y)
    
    Tmatrix <- fields.mkpoly(x, 2)
    qr.T <- qr(Tmatrix)
    N <- nrow(y)
    Q2 <- qr.yq2(qr.T, diag(1, N))
    ys <- t(Q2) %*% y
    
    N2 <- length(ys)
    theta <- exp(ltheta)
    K <- Matern(rdist(x, x)/theta, smoothness = smoothness)
    Ke <- eigen(t(Q2) %*% K %*% Q2, symmetric = TRUE)
    u2 <- t(Ke$vectors) %*% ys
    # mean over replicates   --  mean square for coefficients for
    #  a particular eigenfunction.
    u2.MS <- c(rowMeans(u2^2))
    
    D2 <- Ke$values
    N2 <- length(D2)
    
    
    # grid of lambda based on spacing of eigenvalues
    ngrid <- min(ngrid, N2)
    lambda.grid <- exp(seq(log(D2[1]), log(D2[N2]), , ngrid))
    trA <- minus.pflike <- rep(NA, ngrid)
    #grid search followed by golden section
    
    # -log likelihood
    temp.fn <- function(llam, info) {
        lam.temp <- exp(llam)
        u2 <- info$u2.MS
        D2 <- info$D2
        N2 <- length(u2.MS)
        # MLE of rho
        rho.MLE <- (sum((u2.MS)/(lam.temp + D2)))/N2
        # ln determinant
        lnDetCov <- sum(log(lam.temp + D2))
        -1 * M * (-N2/2 - log(2 * pi) * (N2/2) - (N2/2) * log(rho.MLE) - 
            (1/2) * lnDetCov)
    }
    
    # information list for calling golden section search.
    info <- list(D2 = D2, u2 = u2.MS, M = M)
    out <- golden.section.search(f = temp.fn, f.extra = info, 
        gridx = log(lambda.grid), tol = 1e-07)
    
    minus.LogProfileLike <- out$fmin
    lambda.MLE <- exp(out$x)
    rho.MLE <- (sum((u2.MS)/(lambda.MLE + D2)))/N2
    sigma.MLE <- sqrt(lambda.MLE * rho.MLE)
    trA <- sum(D2/(lambda.MLE + D2))
    pars <- c(rho.MLE, theta, sigma.MLE, trA)
    names(pars) <- c("rho", "theta", "sigma", "trA")
    if (value) {
        return(minus.LogProfileLike)
    }
    else {
        return(list(minus.lPlike = minus.LogProfileLike, lambda.MLE = lambda.MLE, 
            pars = pars, mle.grid = out$coarse.search))
    }
    
}


MLE.Matern.fast <- function(x, y, smoothness, theta.grid = NULL, 
    ngrid = 20, verbose = FALSE, m = 2, ...) {
    # remove missing values and print out a warning
    bad <- is.na(y)
    if (sum(bad) > 0) {
        cat("removed ", sum(bad), " NAs", fill = TRUE)
        x <- x[!bad, ]
        y <- y[!bad]
    }
    
    # list to pass to the objective function
    # NOTE: large ngrid here is very cheap after the eigen decomposition
    #       has been done.
    info <- list(x = x, y = y, smoothness = smoothness, ngrid = 80)
    
    # if grid for ranges is missing use some quantiles of
    #  pairwise distances among data.
    if (is.null(theta.grid)) {
        theta.range <- quantile(rdist(x, x), c(0.03, 0.97))
        theta.grid <- seq(theta.range[1], theta.range[2], , ngrid)
    }
    if (length(theta.grid) == 2) {
        theta.grid <- seq(theta.grid[1], theta.grid[2], , ngrid)
    }
    else {
        ngrid <- length(theta.grid)
    }
    
    # grid search/golden section search
    # note that search is in log scale.
    out <- golden.section.search(f = MLE.objective.fn, f.extra = info, 
        gridx = log(theta.grid))
    
    theta.MLE <- exp(out$x)
    REML <- -out$fmin
    # one final call with the theta.MLE value to recover MLEs for rho and sigma
    out2 <- MLE.objective.fn(log(theta.MLE), info, value = FALSE)
    return(list(smoothness = smoothness, pars = out2$pars[1:3], 
        REML = REML, trA = out2$pars[4], REML.grid = cbind(theta.grid, 
            -1 * out$coarse.search[, 2])))
    
    
}