R/helper.R

In autoFRK: Automatic Fixed Rank Kriging

as.matrix.mrts <- function(x, ...) {
  attr(x, "S") <- NULL
  attr(x, "UZ") <- NULL
  attr(x, "Xu") <- NULL
  attr(x, "nconst") <- NULL
  attr(x, "BBBH") <- NULL
  attr(x, "class") <- NULL
  class(x) <- "matrix"
  return(x)
}

getEigen <- function(A) {
  obj <- getASCeigens(A)
  obj$value <- rev(obj$value)
  obj$vector <- obj$vector[, NCOL(A):1]
  obj
}

getHalf <- function(Fk, iDFk) {
  dec <- getEigen(t(Fk) %*% iDFk)
  dec$vector <- dec$vector
  sroot <- sqrt(pmax(dec$value, 0))
  sroot[sroot == 0] <- Inf
  sroot <- 1 / sroot
  dec$vector %*% (sroot * t(dec$vector))
}

getLikelihood <- function(Data, Fk, M, s, Depsilon) {
  logdet <- function(R, L, K) {
    spam::determinant(diag(1, K) + t(L) %*% solve(R) %*%
                        L, logarithm = TRUE)$modulus + spam::determinant(R,
                                                                         logarithm = TRUE
                        )$modulus
  }
  Data <- as.matrix(Data)
  O <- as.matrix(!is.na(Data))
  TT <- NCOL(Data)
  n2loglik <- sum(O) * log(2 * pi)
  R <- toSpMat(s * Depsilon)
  eg <- eigen(M)
  L <- Fk %*% eg$vector %*% diag(sqrt(pmax(eg$value, 0))) %*%
    t(eg$vector)
  K <- NCOL(Fk)
  for (tt in 1:TT) {
    zt <- Data[O[, tt], tt]
    Rt <- R[O[, tt], O[, tt]]
    Lt <- L[O[, tt], ]
    if (NCOL(Lt) == 1) {
      Lt <- t(Lt)
      n2loglik <- n2loglik + log(Rt + Lt %*% t(Lt))
    }
    else {
      n2loglik <- n2loglik + logdet(Rt, Lt, K) + sum(zt *
                                                       invCz(Rt, Lt, zt))
    }
  }
  return(n2loglik)
}

ifElse <- function(cond, yes_out, no_out) {
  if (cond) {
    return(yes_out)
  } else {
    return(no_out)
  }
}

checkDiag <- function(X) {
  if (is(X, "numeric") & (length(X) == 1)) {
    return(TRUE)
  }
  if (is(X, "matrix")) {
    if (sum(abs(diag(diag(X)) - X)) < .Machine$double.eps) {
      return(TRUE)
    } else {
      return(FALSE)
    }
  }
  else {
    x <- diag.spam(diag.of.spam(X), NROW(X))
    return(identical(x, X))
  }
}

setNC <- function(z, loc, nlevel) {
  Dimension <- NCOL(loc)
  N <- nrow(z)
  a <- sum(2^(Dimension * (0:(nlevel - 1))))
  NCtest <- (N / a)^(1 / Dimension)
  return(round(max(4, NCtest)))
}

subKnot <- function(x, nknot, xrng = NULL, nsamp = 1) {
  x <- as.matrix(x)
  xdim <- dim(x)
  if (xdim[2] > 1) {
    for (kk in xdim[2]:1) x <- x[order(x[, kk]), ]
  }
  else {
    x <- as.matrix(sort(x))
  }
  if (is.null(xrng)) {
    if (xdim[2] > 1) {
      xrng <- apply(x, 2, range)
    }
    else {
      xrng <- matrix(range(x), 2, 1)
    }
  }
  mysamp <- function(zANDid) {
    z <- as.double(names(zANDid))
    if (length(z) == 1L) {
      z
    }
    else {
      set.seed(mean(zANDid))
      sample(z, size = min(nsamp, length(z)))
    }
  }
  rng <- sqrt(xrng[2, ] - xrng[1, ])
  rng[rng == 0] <- min(rng[rng > 0]) / 5
  rng <- rng * 10 / min(rng)
  nmbin <- round(exp(log(rng) * log(nknot) / sum(log(rng))))
  nmbin <- pmax(2, nmbin)
  while (prod(nmbin) < nknot) {
    nmbin[which.max(rng)] <- nmbin[which.max(rng)] +
      1
  }
  gvec <- matrix(1, xdim[1], 1)
  cnt <- 0
  while (length(unique(gvec)) < nknot) {
    nmbin <- nmbin + cnt
    gvec <- matrix(1, xdim[1], 1)
    kconst <- 1
    for (kk in 1:xdim[2]) {
      grp <- pmin(round((nmbin[kk] - 1) * ((x[, kk] - xrng[
        1,
        kk
      ]) / (xrng[2, kk] - xrng[1, kk]))), nmbin[kk] -
        1L)
      if (length(unique(grp)) < nmbin[kk]) {
        brk <- quantile(x[, kk], seq(0, 1, l = nmbin[kk] +
                                       1))
        brk[1] <- brk[1] - 0.1^8
        grp <- as.double(cut(x[, kk], brk))
      }
      gvec <- gvec + kconst * grp
      kconst <- kconst * nmbin[kk]
    }
    cnt <- cnt + 1
  }
  gvec <- as.factor(gvec)
  gid <- as.double(as.character(gvec))
  names(gid) <- 1:xdim[1]
  index <- unlist(tapply(gid, gvec, mysamp))
  if (xdim[2] > 1) {
    x[index, ]
  } else {
    x[index]
  }
}

systemRam <- function(os) {
  remove_white <- function(x) {
    gsub(
      "(^[[:space:]]+|[[:space:]]+$)",
      "", x
    )
  }
  toBytes <- function(value) {
    num <- as.numeric(value[1])
    units <- value[2]
    power <- match(units, c("kB", "MB", "GB", "TB"))
    if (!is.na(power)) {
      return(num * 1024^power)
    }
    power <- match(units, c(
      "Kilobytes", "Megabytes", "Gigabytes",
      "Terabytes"
    ))
    if (!is.na(power)) {
      return(num * 1024^power)
    }
    num
  }
  if (length(grep("^linux", os))) {
    cmd <- "awk '/MemTotal/ {print $2}' /proc/meminfo"
    ram <- system(cmd, intern = TRUE)
    ram <- as.numeric(ram) * 1024
  }
  else if (length(grep("^darwin", os))) {
    ram <- system("system_profiler -detailLevel mini | grep \"  Memory:\"",
                  intern = TRUE
    )[1]
    ram <- remove_white(ram)
    ram <- toBytes(unlist(strsplit(ram, " "))[2:3])
  }
  else if (length(grep("^solaris", os))) {
    cmd <- "prtconf | grep Memory"
    ram <- system(cmd, intern = TRUE)
    ram <- remove_white(ram)
    ram <- toBytes(unlist(strsplit(ram, " "))[3:4])
  }
  else {
    ram <- system("wmic MemoryChip get Capacity", intern = TRUE)[-1]
    ram <- remove_white(ram)
    ram <- ram[nchar(ram) > 0]
    sum(as.numeric(ram))
  }
  as.double(ram)
}

toSpMat <- function(mat) {
  if (is(mat, "data.frame")) {
    mat <- as.matrix(mat)
  }
  if (is(mat, "matrix")) {
    if (length(mat) > 10^8) {
      warnings("Use sparse matrix as input instead; otherwise it could take a very long time!")
      db <- tempfile()
      NR <- NROW(mat)
      NC <- NCOL(mat)
      f <- fm.create(db, NR, NC)
      f[, 1:NCOL(mat)] <- mat
      j <- sapply(1:NC, function(j) which(f[, j] != 0))
      ridx <- unlist(j)
      k <- sapply(1:NR, function(k) rbind(k, which(f[k, ] != 0)))
      kk <- matrix(unlist(k), ncol = 2, byrow = T)
      cidx <- sort(kk[, 2])
      where <- (cidx - 1) * NR + ridx
      closeAndDeleteFiles(f)
    }
    else {
      where <- which(mat != 0)
      ridx <- row(mat)[where]
      cidx <- col(mat)[where]
    }
    nonzero <- mat[where]
    mat <- spam(0, nrow = NROW(mat), ncol = NCOL(mat))
    mat[ridx, cidx] <- nonzero
  }
  return(mat)
}

uniquecombs <- function(x) {
  unique(x)
}

ZinvC <- function(R, L, z) {
  K <- NCOL(L)
  iR <- solve(R)
  ZiR <- z %*% iR
  left <- ZiR %*% L %*% solve(diag(1, K) + t(L) %*% iR %*% L) %*%
    t(L)
  ZiR - left %*% iR
}

print.FRK <- function(x, ...) {
  attr(x, "pinfo") <- NULL
  if (!is.null(x$LKobj)) {
    x$LKobj <- x$LKobj$summary
  }
  out <- paste("a ", NROW(x$G), " by ", NCOL(x$G), " mrts matrix",
               sep = ""
  )
  print(out)
}

print.mrts <- function(x, ...) {
  if (NCOL(x) == 1) {
    out <- c(x)
  } else {
    out <- x[, 1:NCOL(x)]
  }
  print(out)
}

mkpd <- function(M) {
  v <- try(min(eigen(M, only.values = T)$value), silent = TRUE)
  if (is(v, "try-error")) {
    M <- (M + t(M)) / 2
    v <- min(eigen(M, only.values = T)$value)
  }
  if (v <= 0) {
    M <- M + diag(max(0, -v) + 0.1^7.5, NROW(M))
  }
  return(M)
}


extractLK <- function(obj, loc = NULL, w = NULL, pick = NULL) {
  out <- list()
  if (is.null(loc)) {
    loc <- ifElse(
      is.null(obj$LKinfo.MLE$x), obj$LKinfo.MLE$call["x"][[1]],
      obj$LKinfo.MLE$x
    )
  }
  phi <- LKrig.basis(loc, obj$LKinfo)
  Q <- LKrig.precision(obj$LKinfo)
  out$Q <- Q
  if (!is.null(w)) {
    out$weights <- w
  } else {
    out$weights <- obj$LKinfo.MLE$weights
  }
  w <- diag.spam(sqrt(out$weights))
  wX <- w %*% phi
  out$wX <- wX
  out$G <- t(wX) %*% wX + obj$lambda.MLE * Q
  out$lambda <- obj$lambda.MLE
  if (is.null(pick)) {
    pick <- 1:NROW(loc)
  }
  out$pick <- pick
  
  return(out)
}