R/contourutil.R

In excursions: Excursion Sets and Contour Credibility Regions for Random Fields

## Calculate the contour map function
contourfunction.mc <- function(lp, mu, X, ind, alpha, verbose = FALSE) {
  if (missing(mu)) {
    stop("Must supply mu")
  }

  if (missing(lp)) {
    stop("Must supply level plot")
  }

  F.limit <- 1

  if (verbose) cat("set limits\n")

  lim <- excursions.limits(lp = lp, mu = mu, measure = 0)

  m.size <- length(mu)
  indices <- NULL

  if (!missing(ind)) {
    if (is.logical(ind)) {
      indices <- ind
      m.size <- sum(ind)
    } else {
      indices <- rep(FALSE, length(mu))
      indices[ind] <- TRUE
      m.size <- length(ind)
    }
  }
  if (verbose) cat("calculate marginals\n")

  rho <- contourmap.marginals.mc(X = X, lim = lim, ind = indices)


  reo <- excursions.permutation(rho, indices, use.camd = FALSE)

  if (verbose) cat("integrate\n")

  res <- mcint(X = X[reo, ], a = lim$a[reo], b = lim$b[reo])

  n <- length(mu)
  ii <- which(res$Pv[1:n] > 0)
  if (length(ii) == 0) i <- n + 1 else i <- min(ii)

  F <- Fe <- E <- rep(0, n)
  F[reo] <- res$Pv
  Fe[reo] <- res$Ev

  ireo <- NULL
  ireo[reo] <- 1:n

  ind.lowF <- F < 1 - F.limit
  E[F > 1 - alpha] <- 1
  F[ind.lowF] <- Fe[ind.lowF] <- NA

  M <- rep(-1, n)
  for (i in 1:(lp$n.levels + 1)) {
    M[(lp$G == (i - 1)) & (E == 1)] <- i - 1
  }

  return(list(F = F, Fe = Fe, E = E, M = M, rho = rho))
}


## Calculate the contour map function
contourfunction <- function(lp, mu, Q, vars, ind, alpha, n.iter = 10000,
                            F.limit, Q.chol, max.threads = 0,
                            seed = seed, verbose = FALSE,
                            rho, qc = FALSE) {
  if (qc && missing(rho)) {
    stop("Must supply rho if QC method is used.")
  }
  if (!missing(Q.chol) && !is.null(Q.chol)) {
    Q <- Q.chol
    is.chol <- TRUE
  } else if (!missing(Q) && !is.null(Q)) {
    is.chol <- FALSE
  } else {
    stop("Must supply Q or Q.chol")
  }

  if (missing(mu)) {
    stop("Must supply mu")
  }

  if (missing(lp)) {
    stop("Must supply level plot")
  }

  if (missing(F.limit) || is.null(F.limit)) {
    F.limit <- 1
  }

  if (!missing(alpha) && !is.null(alpha)) {
    F.limit <- max(alpha, F.limit)
  }

  if (verbose) cat("set limits\n")
  lim <- excursions.limits(lp = lp, mu = mu, measure = 0)


  if (missing(vars)) {
    if (verbose) cat("calculate variances\n")
    if (is.chol) {
      vars <- excursions.variances(L = Q, max.threads = max.threads)
    } else {
      vars <- excursions.variances(Q = Q, max.threads = max.threads)
    }
  }

  m.size <- length(mu)
  indices <- NULL

  if (!missing(ind)) {
    if (is.logical(ind)) {
      indices <- ind
      m.size <- sum(ind)
    } else {
      indices <- rep(FALSE, length(mu))
      indices[ind] <- TRUE
      m.size <- length(ind)
    }
  }
  if (verbose) cat("calculate marginals\n")
  if (missing(rho) || is.null(rho)) {
    rho <- contourmap.marginals(mu = mu, vars = vars, lim = lim, ind = indices)
    lim$a <- lim$a - mu
    lim$b <- lim$b - mu
  }
  if (qc) {
    lim$a[indices] <- sqrt(vars[indices]) * qnorm(pmin(pmax(rho[indices, 1], 0), 1))
    lim$b[indices] <- sqrt(vars[indices]) * qnorm(pmin(pmax(rho[indices, 2], 0), 1))
    rho <- rho[, 2] - rho[, 1]
  }


  if (verbose) cat("calculate permutation\n")
  use.camd <- !missing(ind) || alpha < 1
  reo <- excursions.permutation(
    rho = rho, ind = indices,
    use.camd = use.camd, alpha = F.limit, Q = Q
  )

  if (verbose) cat("integrate\n")
  res <- excursions.call(lim$a, lim$b, reo, Q,
    is.chol = is.chol,
    1 - F.limit, K = n.iter, max.size = m.size,
    n.threads = max.threads, seed = seed
  )

  n <- length(mu)
  ii <- which(res$Pv[1:n] > 0)
  if (length(ii) == 0) i <- n + 1 else i <- min(ii)

  F <- Fe <- E <- rep(0, n)
  F[reo] <- res$Pv
  Fe[reo] <- res$Ev

  ireo <- NULL
  ireo[reo] <- 1:n

  ind.lowF <- F < 1 - F.limit
  E[F > 1 - alpha] <- 1
  F[ind.lowF] <- Fe[ind.lowF] <- NA

  M <- rep(-1, n)
  for (i in 1:(lp$n.levels + 1)) {
    M[(lp$G == (i - 1)) & (E == 1)] <- i - 1
  }

  return(list(F = F, Fe = Fe, E = E, M = M, rho = rho))
}

## Calculate marginal probabilities P(lim$a < X < lim$b) for
## when X is N(mu,vars)-distributed
contourmap.marginals <- function(mu, vars, lim, ind) {
  if (!missing(ind) && !is.null(ind)) {
    marg <- rep(0, length(mu))
    marg[ind] <- pnorm(lim$b[ind], mu[ind], sqrt(vars[ind])) -
      pnorm(lim$a[ind], mu[ind], sqrt(vars[ind]))
  } else {
    marg <- pnorm(lim$b, mu, sqrt(vars)) - pnorm(lim$a, mu, sqrt(vars))
  }
  return(marg)
}

## Calculate marginal probabilities P(lim$a < X < lim$b) for
## when X
contourmap.marginals.mc <- function(X, lim, ind) {
  if (!missing(ind) && !is.null(ind)) {
    marg <- rep(0, length(lim$a))
    marg[ind] <- rowMeans(lim$a[ind] < X[ind, ] & X[ind, ] < lim$b[ind])
  } else {
    marg <- rowMeans(lim$a < X & X < lim$b)
  }
  return(marg)
}

## Create a levelplot with given levels/number of levels
excursions.levelplot <- function(mu, n.levels, ind, levels,
                                 equal.area = FALSE,
                                 pretty.cm = FALSE) {
  n <- length(mu)
  if (missing(ind)) ind <- 1:n
  x.mean <- rep(-Inf, n)
  x.mean[ind] <- mu[ind]
  r <- range(mu[ind])
  if (missing(levels)) {
    if (missing(n.levels)) {
      stop("Must specify levels or number of levels")
    } else {
      private.check.integer(n.levels)
    }
    if (equal.area) {
      levels <- as.vector(quantile(mu[ind],
        (1:n.levels) / (n.levels + 1),
        type = 4
      ))
    } else if (pretty.cm) {
      levels <- pretty(mu[ind], n = n.levels)
      n.levels <- length(levels)
    } else {
      levels <- seq(
        from = r[1], to = r[2],
        length.out = (n.levels + 2)
      )[2:(n.levels + 1)]
    }
  } else {
    if (!missing(n.levels)) {
      if (n.levels != length(levels)) {
        warning("n.levels is not equal to the length of levels")
        n.levels <- length(levels)
      }
    } else {
      n.levels <- length(levels)
    }
  }

  u.e <- NULL
  E <- vector("list", n.levels + 1)

  if (n.levels == 1) {
    l1 <- c(levels - diff(r), levels, levels + diff(r))
  } else {
    l1 <- c(
      2 * levels[1] - levels[2],
      levels,
      2 * levels[n.levels] - levels[n.levels - 1]
    )
  }
  if (min(x.mean[ind]) < l1[1]) {
    l1[1] <- min(min(x.mean[x.mean > -Inf]), l1[1])
  }

  for (i in 1:(n.levels + 1)) u.e[i] <- (l1[i] + l1[i + 1]) / 2

  for (i in 1:(n.levels)) {
    E[[i]] <- which((l1[i] <= x.mean) & (x.mean < l1[i + 1]))
  }
  E[[n.levels + 1]] <- which((l1[n.levels + 1] <= x.mean))
  map <- rep(0, n)
  G <- rep(-1, n)
  for (i in 1:(n.levels + 1)) {
    map[E[[i]]] <- u.e[i]
    G[E[[i]]] <- i - 1
  }
  return(list(
    u = levels,
    n.levels = n.levels,
    u.e = u.e,
    E = E,
    map = map,
    G = G
  ))
}

## Create a P-optimal levelplot.
## The function will take A LOT of time to run if use.marginals=FALSE.
excursions.opt.levelplot <- function(mu, vars, Q, n.levels, measure = 2, 
                                     use.marginals = TRUE, ind, max.threads = 0) {
  if ((measure != 1) && (measure != 2) && (measure != 0)) {
    stop("only measure 0, 1, or 2 allowed")
  }

  if (missing(ind)) ind <- seq_along(mu)

  r <- range(mu[ind])

  u <- seq(from = r[1], to = r[2], length.out = (n.levels + 2))[2:(n.levels + 1)]
  if (n.levels == 1) {
    l <- diff(r)
  } else {
    l <- diff(u[1:2])
  }

  Q.chol <- chol(Q)
  u.add <- seq(from = -l, to = l, length.out = 19)
  P.add <- NULL
  for (jj in seq_along(u.add)) {
    P.add[jj] <- -restricted.lim.func(
      u.add = u.add[jj],
      u0 = u,
      mu = mu,
      vars = vars,
      Q.chol = Q.chol,
      Q = Q,
      measure = measure,
      use.marginals = use.marginals,
      ind = ind,
      max.threads = max.threads
    )
  }
  plot(u.add, P.add)
  k <- which.max(P.add)
  u.add <- u.add[k]
  P.k <- P.add[k]

  u <- u.add + u
  lp <- excursions.levelplot(mu, levels = u, ind = ind)
  if (measure == 2) {
    if (use.marginals) {
      lp$P2.bound <- P.k
    } else {
      lp$P2 <- P.k
    }
  } else if (measure == 1) {
    if (use.marginals) {
      lp$P1.bound <- P.k
    } else {
      lp$P1 <- P.k
    }
  } else if (measure == 0) {
    if (use.marginals) {
      lp$P0.bound <- P.k
    } else {
      lp$P0 <- P.k
    }
  }
  return(lp)
}
## Internal function for optimization of restricted P-optimal contour map
restricted.lim.func <- function(u.add, u0, mu, vars, Q.chol, Q, measure,
                                use.marginals, ind = ind, max.threads = 0) {
  lp <- excursions.levelplot(mu = mu, levels = (u.add + u0), ind = ind)
  if (use.marginals) {
    val <- -Pmeasure.bound(lp = lp, mu = mu, vars = vars, type = measure, ind = ind)
    cat(u0, ": ", -val, "\n")
  } else {
    val <- -Pmeasure(lp,
      mu = mu, Q = Q, Q.chol = Q.chol, type = measure,
      ind = ind, vars = vars, max.threads = max.threads
    )
    cat(u.add, ": ", -val, "\n")
  }
  return(val)
}


## Internal function for optimization of P-optimal contour map
excursions.lim.func <- function(u, mu, vars, Q.chol, Q, measure,
                                use.marginals, ind = ind, max.threads = 0) {
  lp <- excursions.levelplot(mu, levels = u, ind = ind)
  if ((min(u) <= min(mu[ind])) || (max(u) >= max(mu[ind]))) {
    # levels should be in (min(mu),max(mu))
    val <- 0
  } else if (max(sort(u, index.return = TRUE)$ix - seq_len(length(u))) > 0) {
    # levels should be sorted
    val <- 0
  } else {
    v <- TRUE
    if (length(u) > 1) {
      for (i in 1:(length(u) - 1)) {
        v <- v & (sum(mu[ind] < u[i + 1] & mu[ind] > u[i]) > 0)
      }
    }
    if (!v) {
      # all sets E should be non-empty
      val <- 0
    } else {
      if (use.marginals) {
        val <- -Pmeasure.bound(lp = lp, mu = mu, vars = vars, type = measure, ind = ind)
        cat(u, ": ", -val, "\n")
      } else {
        val <- -Pmeasure(lp,
          mu = mu, Q = Q, Q.chol = Q.chol, type = measure,
          ind = ind, vars = vars, max.threads = max.threads
        )
        cat(u, ": ", -val, "\n")
      }
    }
  }
  return(val)
}

Pmeasure.bound <- function(lp, mu, vars, type, ind = NULL) {
  limits <- excursions.limits(lp, mu, measure = type)
  if (type == 0) {
    return(mean(contourmap.marginals(mu, vars, limits, ind)[ind]))
  } else {
    return(min(contourmap.marginals(mu, vars, limits, ind)[ind]))
  }
}

## Function that calculates the P measure for a given contour map.
Pmeasure <- function(lp, mu, Q, Q.chol, ind = NULL, type, 
                     vars = vars, seed = NULL, n.iter = NULL,
                     max.threads = 0) {
  if (type == 0) {
    res <- contourfunction(lp = lp, mu = mu, Q = Q, vars = vars, ind = ind,
                            max.threads = max.threads)
    p <- mean(res$F[ind])
  } else {
    if (type == 1 && length(lp$u) == 1) {
      return(1)
    }
    limits <- excursions.limits(lp = lp, mu = mu, measure = type)
    res <- gaussint(
      mu = mu, Q = Q, Q.chol = Q.chol, a = limits$a,
      b = limits$b, ind = ind, use.reordering = "limits",
      n.iter = n.iter, seed = seed,
      max.threads = max.threads
    )
    p <- res$P[1]
  }
  return(list(P = res$P[1], E = res$E[1]))
}

## Function that calculates the P measure for a given contour map.
Pmeasure.mc <- function(lp, mu, X, ind = NULL, type) {
  if (type == 0) {
    res <- contourfunction.mc(lp = lp, X = X, ind = ind)
    p <- mean(res$F[ind])
  } else {
    if (type == 1 && length(lp$u) == 1) {
      return(1)
    }
    limits <- excursions.limits(lp = lp, mu = mu, measure = type)
    res <- mcint(X = X, a = limits$a, b = limits$b, ind = ind)
    p <- res$P[1]
  }
  return(p)
}

## Set integration limits for a given measure.
excursions.limits <- function(lp, mu, measure) {
  n <- length(mu)
  n.l <- length(lp$u)

  a <- rep(-Inf, n)
  b <- rep(Inf, n)

  if (measure == 2 || measure == "P2" || measure == "P2-bound") {
    b[lp$E[[1]]] <- lp$u.e[2]
    a[lp$E[[n.l + 1]]] <- lp$u.e[n.l]

    if (n.l > 1) {
      for (i in 2:n.l) {
        a[lp$E[[i]]] <- lp$u.e[i - 1]
        b[lp$E[[i]]] <- lp$u.e[i + 1]
      }
    }
  } else if (measure == 1 || measure == "P1" || measure == "P1-bound") {
    if (n.l > 1) {
      b[lp$E[[1]]] <- lp$u[2]
      a[lp$E[[n.l + 1]]] <- lp$u[n.l - 1]
    }
    if (n.l > 2) {
      b[lp$E[[2]]] <- lp$u[3]
      a[lp$E[[n.l]]] <- lp$u[n.l - 2]
    }
    if (n.l > 3) {
      for (i in 3:(n.l - 1)) {
        a[lp$E[[i]]] <- lp$u[i - 2]
        b[lp$E[[i]]] <- lp$u[i + 1]
      }
    }
  } else if (measure == 0 || measure == "P0" || measure == "P0-bound") {
    b[lp$E[[1]]] <- lp$u[1]
    a[lp$E[[n.l + 1]]] <- lp$u[n.l]

    if (n.l > 1) {
      for (i in 2:n.l) {
        a[lp$E[[i]]] <- lp$u[i - 1]
        b[lp$E[[i]]] <- lp$u[i]
      }
    }
  } else {
    stop("Measure must be 0, 1, or 2")
  }
  return(list(a = a, b = b))
}

#' Define a color map for displaying contour maps.
#'
#' \code{contourmap.colors} calculates suitable colours for displaying contour maps.
#'
#' @param lp A contourmap calculated by \code{contourmap}, \code{contourmap.inla}, or \code{contourmap.mc}
#' @param zlim The range that should be used (optional). The default is the range of the mean value function used when creating the contourmap.
#' @param col The colormap that the colours should be taken from.
#' @param credible.col The color that should be used for displaying the credible regions for the contour curves (optional).
#'
#' @return A color map.
#' @author David Bolin \email{davidbolin@@gmail.com}
#' @export
#'
#' @examples
#' n <- 10
#' Q <- Matrix(toeplitz(c(1, -0.5, rep(0, n - 2))))
#' map <- contourmap(
#'   mu = seq(-5, 5, length = n), Q, n.levels = 2,
#'   compute = list(F = FALSE), max.threads = 1
#' )
#' cols <- contourmap.colors(map,
#'   col = heat.colors(100, 1),
#'   credible.col = grey(0.5, 1)
#' )
contourmap.colors <- function(lp, zlim, col, credible.col) {
  if (missing(zlim) || is.null(zlim)) {
    zlim <- lp$meta$mu.range
  }

  u.e <- sort(unique(lp$map)) # extracts only used breakpoints
  u.e[1] <- max(u.e[1], zlim[1])
  u.e[length(u.e)] <- min(u.e[length(u.e)], zlim[2])

  breaks <- seq(zlim[1] - 1e-16, zlim[2] + 1e-16,
    length.out = length(col) + 1
  )
  cmap <- col[cut(c(u.e), breaks)@.Data]
  if (!missing(credible.col) && !is.null(credible.col)) {
    cmap <- c(credible.col, cmap)
  }

  return(cmap)
}