#### create SGP emulator ####

## Do not edit this file, which is created by 'notangle SGP.nw > SGP.R'.
## Edit SGP.nw instead.


#### set up emulator using prior hyperparameters and ensemble

"setUpSGP" <- function(g, r, m, V, a, d, f = NULL, X = NULL, ir = NULL)
{
  verbose <- getOption('verbose')
  verbose <- is.logical(verbose) && verbose

  ## unpack list argument

  if (inherits(g, "SGP")) {
    SGP <- g
    for(obj in c("g", "r", "m", "V", "a", "d", "ir"))
      assign(obj, SGP[[obj]])
    rm(SGP)
  }

  ## package into a list: initial objects

  robj <- list(g = g, r = r, m = m, V = V, a = a, d = d,
    f = f, X = X, ir = ir)

  ## check the arguments

  stopifnot(is.vector(m), is.matrix(V), length(a) == 1, length(d) == 1)
  q <- length(m)
  stopifnot(identical(dim(V), c(q, q)))

  ## chance for a quick exit

  if (is.null(f) || is.null(X)) {

    class(robj) <- c("SGP", class(robj))
    return(robj)

  } else {

    stopifnot(is.function(g), is.function(r))

    G <- try(g(X))
    if (inherits(G, "try-error"))
      stop("Cannot evaluate 'g' with argument 'X'")
    stopifnot(is.matrix(G), ncol(G) == q, !any(is.na(G)))

    n <- nrow(G)
    stopifnot(is.vector(f), length(f) == n)

    Rx <- try(r(X, X))
    if (inherits(Rx, "try-error"))
      stop("Cannot evaluate 'r' with argument 'X, X'")
    stopifnot(is.matrix(Rx), identical(dim(Rx), c(n, n)))

  }

  if (is.null(ir)) {

    ## update the B mean and scale

    Q <- chol(crossprod(tcrossprod(chol(V), G)) + Rx, pivot = TRUE) # chol(S)
    pivot <- attr(Q, "pivot")

    PP <- backsolve(Q, cbind(G %*% V, f - G %*% m)[pivot, , drop = FALSE],
      transpose = TRUE)
    P1 <- PP[, 1:q, drop = FALSE]
    P2 <- PP[, q + 1]
    rm(PP)

    EB <- m + drop(crossprod(P1, P2))
    SB <- V - crossprod(P1)

    ## update a and d

    anew <- a + n
    dnew <- d + drop(crossprod(P2))

    ## add new objects

    robj <- c(robj, list(Q = Q, P1 = P1, P2 = P2))

  } else {

    ## handle non-null ir

    if (!is.null(ir)) {

      if (!is.function(ir)) {
        if (verbose)
          warning("Providing default 'ir' function")
        "ir" <- function(X) chol2inv(chol(r(X, X)))
      }
    }

    ## update hyperparameters using SMW approach

    iRx <- try(ir(X))
    if (inherits(iRx, "try-error"))
      stop("Cannot evaluate 'ir' with argument 'X'")
    stopifnot(is.matrix(iRx), identical(dim(iRx), c(n, n)))

    Q <- chol(chol2inv(chol(V)) + crossprod(G, iRx %*% G)) # no decomp of iRx
    SB <- chol2inv(Q)

    iS <- iRx - crossprod(backsolve(Q, crossprod(G, iRx), transpose = TRUE))
    P1 <- tcrossprod(V, G)
    P2 <- drop(iS %*% (f - G %*% m))

    EB <- m + drop(P1 %*% P2)

    anew <- a + n
    dnew <- d + drop(crossprod(f - G %*% m, P2))

    ## add new objects

    robj$ir <- ir # might have been redefined
    robj <- c(robj, list(iS = iS, P1 = P1, P2 = P2))
    
  }

  ## add updated hyperparameters and return

  robj <- c(robj, list(EB = EB, SB = SB, anew = anew, dnew = dnew))
  class(robj) <- c("SGP", class(robj))
  return(robj)
}


"predictSGP" <- function(SGP, Z, type = c("pars", "MV"),
  fromPrior = FALSE, r2 = r)
{
  stopifnot(inherits(SGP, "SGP"), is.logical(fromPrior))
  type <- match.arg(type)
  verbose <- getOption('verbose')
  verbose <- is.logical(verbose) && verbose

  ## unpack basic objects

  for (obj in c("g", "r", "f", "X"))
    assign(obj, SGP[[obj]])

  if (!identical(r, r2))
    r <- r2

  ## option to predict with the prior

  if (fromPrior || is.null(f) || is.null(X)) {

    if (verbose)
      warning("Predicting from the prior")

    for (obj in c("m", "V", "a", "d"))
      assign(obj, SGP[[obj]])

    X <- Z
    f <- rep(0, nrow(X))
    q <- length(m)

    stopifnot(is.function(g), is.function(r))

    G <- try(g(X))
    if (inherits(G, "try-error"))
      stop("Cannot evaluate 'g' with argument 'X'")
    stopifnot(is.matrix(G), ncol(G) == q, !any(is.na(G)))

    n <- nrow(G)
    stopifnot(is.vector(f), length(f) == n)

    Rx <- try(r(X, X))
    if (inherits(Rx, "try-error"))
      stop("Cannot evaluate 'r' with argument 'X, X'")
    stopifnot(is.matrix(Rx), identical(dim(Rx), c(n, n)))

    ## do prediction with prior

    EfZ <- drop(G %*% m)
    SfZ <- crossprod(tcrossprod(chol(V), G)) + Rx
    anew <- a
    dnew <- d

    ## package and return

    if (type == "pars")
      return(list(E = EfZ, S = (dnew / anew) * SfZ, df = anew))
    else if (type == "MV")
      return(list(M = EfZ, V = (dnew / (anew - 2)) * SfZ))
    else stop("Never get here.")

  }

  ## predicting from the updated distribution

  if (is.null(SGP$ir)) {

    ## unload new objects

    for (obj in c("Q", "P1", "P2"))
      assign(obj, SGP[[obj]])

    ## compute mean and scale of uz

    Rz <- r(Z, Z)
    Rxz <- r(X, Z)

    pivot <- attr(Q, "pivot")
    P3 <- backsolve(Q, Rxz[pivot, , drop = FALSE], transpose = TRUE)
    EuZ <- drop(crossprod(P3, P2))
    SuZ <- Rz - crossprod(P3)
    SZB <- -crossprod(P3, P1)

  } else {

    ## unload new objects

    for (obj in c("iS", "P1", "P2"))
      assign(obj, SGP[[obj]])

    ## compute mean and variance of uz

    Rz <- r(Z, Z)
    Rxz <- r(X, Z)

    EuZ <- drop(crossprod(Rxz, P2))
    P3 <- iS %*% Rxz
    SuZ <- Rz - crossprod(Rxz, P3) # no decomp of iS
    SZB <- -crossprod(P3, t(P1))

  }

  ## Conditional mean and scale of fz

  for (obj in c("EB", "SB", "anew", "dnew"))
    assign(obj, SGP[[obj]])

  G <- try(g(Z))
  if (inherits(G, "try-error"))
    stop("Problem with evaluating 'g' at 'Z', check 'Z'")

  EfZ <- drop(G %*% EB) + EuZ
  SfZ <- crossprod(tcrossprod(chol(SB), G)) + SuZ
  tmp <- tcrossprod(SZB, G)
  SfZ <- SfZ + tmp + t(tmp)

  ## package and return

  if (type == "pars")
    return(list(E = EfZ, S = (dnew / anew) * SfZ, df = anew))
  else if (type == "MV")
    return(list(M = EfZ, V = (dnew / (anew - 2)) * SfZ))
  else stop("Never get here.")

}

#### Simple example, run with 'eval(SGPexample)'

# source("SGP.R")

"SGPexample" <- expression({
  ## Create a function on [0, 1]

  foo <- function(x) 1 + (3/2) * x + exp(-x) * sin(3 * pi * x / 2)

  ## sample some points

  n <- 5
  x <- runif(n)
  f <- foo(x)

  ## Set up an emulator

  "g" <- function(x) {
    v <- 2*x - 1 # onto [-1, 1]
    cbind(1, v, v^2) # constant, linear, and quadratic terms
  }

  theta <- 0.33 / sqrt(-log(0.05)) # correlation length 0.33
  "r" <- function(x, y) {
    x <- as.vector(x)
    y <- as.vector(y)
    exp(-(outer(x, y, "-")/ theta)^2)
  }

  myEm <- setUpSGP(
    g = g, r = r,
    m = c(1, 0, 0), V = diag(0.1, 3), a = 3, d = 1,
    f = f, X = matrix(x))

  ## make a prediction along a grid

  grid <- seq(from = 0, to = 1, length = 101)
  pp <- predictSGP(myEm, Z = matrix(grid), type = "pars")

  ## draw a picture?

  if (interactive()) {

    ## little function to draw picture

    "drawPicture" <- function(SGP, grid, fromPrior = FALSE, main = "",
      alpha = 0.05)
    {
      pp <- predictSGP(SGP, Z = matrix(grid), type = "pars",
        fromPrior = fromPrior)

      sd <- sqrt(diag(pp$S))
      z <- qt(c(alpha/2, 1 - alpha/2), df = pp$df)
      lims <- pp$E + outer(sd, z); colnames(lims) <- c("lower", "upper")

      mu <- if (fromPrior)
        drop(g(matrix(grid)) %*% SGP$m)
      else
        drop(g(matrix(grid)) %*% SGP$EB)

      plot(foo, type = "n", axes = FALSE, ylim = range(lims, mu),
        main = main, cex.main = 1)

      polygon(c(grid, rev(grid)), c(lims[, "lower"], rev(lims[, "upper"])),
        col = "grey", border = NA)

      axis(1)
      axis(2)

      if (!fromPrior)
        points(as.vector(SGP$X), SGP$f, pch = 16)
      lines(grid, mu, lty = 2)

      invisible(NULL)
    }

    ## from the prior

    op <- par(ask = TRUE)

    drawPicture(myEm, grid = grid, fromPrior = TRUE,
      main = "Function, prior marginal 95% CI and regression line")
    lines(grid, foo(grid), lwd = 2)

    ## from the updated distribution

    drawPicture(myEm, grid = grid, fromPrior = FALSE,
      main = "Function, evaluations, updated marginal 95% CI and regression line")
    lines(grid, foo(grid), lwd = 2)

    par(op)
  }

  ## sample from the emulator

  sam <- rsamSGP(myEm, Z = matrix(grid), n = 10)
  if (interactive()) {
    matplot(grid, t(sam), add = TRUE, type = "l", lty = 1)
    if (!file.exists("SGPpicture.pdf"))
      dev.print(pdf, file = "SGPpicture.pdf")
  }
  ## Compute the negative logarithmic score

  nls <- nlsSGP(myEm)
  cat(sprintf("Negative log score is %.2f\n", nls))

  ## compute subset diagnostics

  loo <- looSGP(myEm)
  cat("Leave-one-out diagnostic:\n")
  print(loo)

  osa <- osaSGP(myEm)
  cat("One-step-ahead diagnostic:\n")
  print(osa)

  stopifnot(all.equal.numeric(nls, sum(osa)))

  #### try out the SMW emulator

  myEm1 <- setUpSGP(
    g = g, r = r,
    m = c(1, 0, 0), V = diag(0.1, 3), a = 3, d = 1,
    f = f, X = matrix(x), ir = NA) # trigger SMW with non-null ir

  pp1 <- predictSGP(myEm1, Z = matrix(grid))

  stopifnot(all.equal(pp, pp1)) # tough test!

  if (interactive()) {
    op <- par(ask = TRUE)
    drawPicture(myEm1, grid = grid,
    main = "SMW Emulator (different updating equations, same result)")
    par(op)
  }

})

#### Sample from the emulator

## MV Student-t random variates

"mvrt" <- function(n = 1, m, S, delta, tol = 1e-6, drop = FALSE)
{
  k <- length(m)
  stopifnot(is.matrix(S), identical(dim(S), c(k, k)), length(delta) == 1)

  eS <- eigen(S, symmetric = TRUE)
  ev <- eS$values
  if (!all(ev >= -tol * abs(ev[1])))
      stop("'S' is not non-negative definite")
  X <- matrix(rnorm(k * n), n)

  X <- eS$vectors %*% diag(sqrt(pmax(ev, 0)), k) %*% t(X)
  chi <- rchisq(n, df = delta)
  X <- drop(m) + sweep(X, 2, sqrt(chi/delta), "/")

  nm <- names(m)
  if (is.null(nm) && !is.null(dn <- dimnames(S)))
      nm <- dn[[1]]
  dimnames(X) <- list(nm, NULL)
  if (n == 1 && drop)
    drop(X)
  else t(X)
}

## wrapper for the emulator

"rsamSGP" <- function(SGP, Z, n = 1, drop = (n == 1), fromPrior = FALSE)
{
  stopifnot(inherits(SGP, "SGP"), length(n) == 1, n >= 1)
  pp <- predictSGP(SGP, Z = Z, type = "pars", fromPrior = fromPrior)
  mvrt(n = n, m = pp$E, S = pp$S, delta = pp$df, drop = drop)
}

#### Emulator diagnostics

## neg log score

"nlsSGP" <- function(SGP, fz = NULL, Z = NULL)
{
  if (is.null(fz) || is.null(Z)) {

    ## check for ensemble

    stopifnot(inherits(SGP, "SGP"))
    if (is.null(SGP$f) || is.null(SGP$X))
      stop("Cannot proceed without an ensemble!")

    fz <- SGP$f
    Z <- SGP$X
    fromPrior <- TRUE

  } else fromPrior <- FALSE

  ## get prior predictive distribution

  pp <- predictSGP(SGP, Z = Z, type = "pars", fromPrior = fromPrior)
  k <- length(pp$E)
  delta <- pp$df

  ## compute negative of [1 + ...]^{...}

  Q <- try(chol(pp$S))
  if (inherits(Q, "try-error")) {
    warning("Not positive definite scale matrix")
    return(NA)
  }
  q <- drop(crossprod(backsolve(Q, fz - pp$E, transpose  = TRUE)))
  robj <- ((delta + k)/2) * log(1 + q / delta)

  ## add neg normalising constant

  robj <- robj + sum(log(diag(Q))) + (k/2) * log(delta * pi) -
    lgamma((delta + k) / 2) + lgamma(delta / 2)
  robj
}

## Leave-one-out

"looSGP" <- function(SGP)
{
  ## check for ensemble

  stopifnot(inherits(SGP, "SGP"))
  if (is.null(SGP$f) || is.null(SGP$X))
    stop("Cannot proceed without an ensemble!")

  loo <- lapply(seq(along = SGP$f), function(i) {
    em <- setUpSGP(SGP, f = SGP$f[-i], X = SGP$X[-i, , drop = FALSE])
    nlsSGP(em, fz = SGP$f[i], Z = SGP$X[i, , drop = FALSE])
  })

  unlist(loo)
}

## cross-validation score

"cvsSGP" <- function(SGP) sum(looSGP(SGP))

## One-step-ahead

"osaSGP" <- function(SGP)
{
 ## check for ensemble

 stopifnot(inherits(SGP, "SGP"))
 if (is.null(SGP$f) || is.null(SGP$X))
   stop("Cannot proceed without an ensemble!")

  n <- nrow(SGP$X)
  osa <- vector("list", n)

  for (i in 1:n) {

    if (i == 1)
      em <- setUpSGP(SGP)
 
    osa[[i]] <- nlsSGP(em, fz = SGP$f[i], Z = SGP$X[i, , drop = FALSE])

    if (i < n)
      em <- setUpSGP(SGP, f = SGP$f[1:i], X = SGP$X[1:i, , drop=FALSE])
  }
 
  unlist(osa)
}


#### fit correlation length by x-validation

"SGPexampleFit" <- expression({
  ## make objective function of correlation length

  "optimand" <- function(theta, type = c("CVS", "NLS")) {

    type <- match.arg(type)

    ## r inherits local theta

    "r" <- function(x, y) {
      x <- as.vector(x)
      y <- as.vector(y)
      exp(-(outer(x, y, "-")/ theta)^2)
    }

    em <- setUpSGP(g = myEm$g, r = r, m = myEm$m, V = myEm$V,
      a = myEm$a, d = myEm$d, f = myEm$f, X = myEm$X)

    if (type == "CVS")
      cvsSGP(em)
    else if (type == "NLS")
      nlsSGP(em)
    else stop("Never get here.")
  }

  ## call optimize

  opt <- optimize(f = optimand,
    lower = 0.1 / sqrt(-log(0.05)),  # correlation length 0.1
    upper = 1.0 / sqrt(-log(0.05)),  # correlation length 1.0
    type = "CVS")

  cat(sprintf("Optimal value is theta = %.2f; corr length = %.2f\n",
    opt$minimum, opt$minimum * sqrt(-log(0.05))))

  ## Note that this solution is likely to be on the upper boundary

  if (interactive()) {
    op <- par(ask = TRUE)
    theta <- opt$minimum
    myEmOpt <- setUpSGP(myEm, f = myEm$f, X = myEm$X) # r() picks up new theta
    drawPicture(myEmOpt, grid = grid, fromPrior = FALSE,
      main = "Function, evaluations, updated marginal 95% CI and regression line\nAfter tuning the correlation length")
    lines(grid, foo(grid), lwd = 2)
    par(op)
  }

  ## Let this be a cautionary tale about the difficulty of interpreting
  ## the regression coefficients when the correlation length is fitted!

})

#### A more complicated example using ir, run with 'eval(SGPexampleIR)'

"SGPexampleIR" <- expression({

  ## spherical stationary isotropic correlation function
  ## with correlation length a (corr is zero beyond a).

  "spherical" <- function(h, a) {
    dd <- dim(h)
    h <- pmin(a, abs(h))
    robj <- 1 - (3/2) * (h / a) + (1/2) * (h / a)^3
    dim(robj) <- dd
    robj
  }

  ## correlation function can handle product structure in X

  clen <- list("locs" = 10, "times" = 1)

  "r" <- function(x, y) {

    if (is.call(x) && x[[1]] == "expand.grid" && identical(y, x)) {

      Vlocs <- spherical(outer(x$locs, x$locs, "-"), a = clen$locs)
      Vtimes <- spherical(outer(x$times, x$times, "-"), a = clen$times)

      return(kronecker(Vtimes, Vlocs))
    }

    if (is.call(x)) x <- eval(x)
    if (is.call(y)) y <- eval(y)

    Vlocs <- spherical(outer(x$locs, y$locs, "-"), a = clen$locs)
    Vtimes <- spherical(outer(x$times, y$times, "-"), a = clen$times)

    return(Vtimes * Vlocs)
  }

  ## Inverse correlation function insists on product structure

  ## Pass in SparseM flag through getOption()

  "ir" <- function(x) {

    useSparseM <- getOption('useSparseM')
    useSparseM <- is.logical(useSparseM) && useSparseM

    if (is.call(x) && x[[1]] == "expand.grid") {

      Vlocs <- spherical(outer(x$locs, x$locs, "-"), a= clen$locs)
      Vtimes <- spherical(outer(x$times, x$times, "-"), a = clen$times)

      if (useSparseM) {

        stopifnot(require(SparseM))

        Vlocs <- as.matrix.csr(Vlocs)   # could be more
        Vtimes <- as.matrix.csr(Vtimes) # sophisticated here!

        return(as.matrix(solve(Vtimes) %x% solve(Vlocs)))

      } else {

        return(chol2inv(chol(Vtimes)) %x% chol2inv(chol(Vlocs)))
      }

    } else stop('Unexpected argument')
  }

  ## make up some data and build an emulator

  X <- call("expand.grid",
         locs = 1:30,
         times = seq(from = 0, to = 5, by = 1/12))

  tmp <- eval(X)
  f <- tmp$locs + 0.5 * (tmp$locs - mean(tmp$locs))^2 +
    tmp$locs * sin(2 * pi * tmp$times)
  rm(tmp)

  "g" <- function(X) {
    X <- eval(X)
    cbind(1, X$locs, X$times, X$locs * X$times, X$locs^2)
  }

  cat("system.time() for different treatments ...\n")

  op <- options(useSparseM = NULL)

  cat("Default, no ir:\n")
  print(system.time(
    myEm2 <- setUpSGP(g = g, r = r, m = rep(0, 5), V = diag(1, 5),
      a = 3, d = 1, f = f, X = X)
  ))

  cat("With default ir:\n")
  print(system.time(
    myEm3 <- setUpSGP(g = g, r = r, m = rep(0, 5), V = diag(1, 5),
      a = 3, d = 1, f = f, X = X, ir = NA)
  ))

  cat("With specified ir:\n")
  print(system.time(
    myEm4 <- setUpSGP(g = g, r = r, m = rep(0, 5), V = diag(1, 5),
      a = 3, d = 1, f = f, X = X, ir = ir)
  ))

  if (require(SparseM)) {
    cat("With specified ir, and SparseM:\n")
    options(useSparseM = TRUE)
    print(system.time(
      myEm5 <- setUpSGP(g = g, r = r, m = rep(0, 5), V = diag(1, 5),
        a = 3, d = 1, f = f, X = X, ir = NA)
    ))
  }

  options(op)

  ## Check they all give the same predictions

  cat("Checking the answers are all the same ...\n")

  Z <- call("expand.grid",
         locs = 1:3 + 0.5,
         times = c(2, 4, 6))

  pp2 <- predictSGP(myEm2, Z = Z)
  pp3 <- predictSGP(myEm3, Z = Z); stopifnot(all.equal(pp2, pp3))
  pp4 <- predictSGP(myEm4, Z = Z); stopifnot(all.equal(pp3, pp4))
  if (exists("myEm5", envir = .GlobalEnv, mode = "list")) {
    pp5 <- predictSGP(myEm5, Z = Z); stopifnot(all.equal(pp4, pp5))
  }

})


#### option to run examples if debug == TRUE

debug <- getOption('debug')
if (is.logical(debug) && debug) {
  eval(SGPexample)
  eval(SGPexampleFit)
  eval(SGPexampleEB)
  eval(SGPexampleMix)
  eval(SGPexampleIR)
}

#### choose non-linear transformation to plug-in, run with
#### 'eval(SGPexampleEB)'

"SGPexampleEB" <- expression({

  ## create a Box-Cox function

  "BoxCox" <- function(x, lambda = 1, tol = 1e-3) {
    if (abs(lambda) < tol)
      log(x)
    else (x^lambda - 1) / lambda
  }

  ## wrapper to transform SGP with a value for lambda

  "BoxCoxSGP" <- function(SGP, lambda)
  {
    stopifnot(inherits(SGP, "SGP"))

    f <- SGP$f
    logJ <- (lambda - 1) * sum(log(f))
    ell <- exp(logJ / length(f))

    f <- BoxCox(f, lambda)

    em <- setUpSGP(g = SGP$g, r = SGP$r, m = SGP$m * ell, V = SGP$V, a = SGP$a,
      d = SGP$d * ell^2, f = f, X = SGP$X, ir = SGP$ir)

    nls <- nlsSGP(em) - logJ
    attr(em, "BoxCox") <- list(lambda = lambda, nls = nls)
    
    em
  }

  ## objective function, starting from emulator Em0

  "nlsLam" <- function(lambda) {
    em <- BoxCoxSGP(Em0, lambda = lambda)
    attr(em, "BoxCox")$nls
  }

  ## minimise on lambda in [-3, 3]

  if (!exists("myEm", inherits = FALSE))
    stop("Do 'eval(SGPexample)' first")

  Em0 <- myEm
  opt <- optimize(f = nlsLam, interval = c(-3, 3))

  ## But more fun to draw the picture

  if (interactive()) {
    op <- par(ask = TRUE)
    lambda <- seq(from = -3, to = 3, by = 0.5)
    y <- sapply(lambda, nlsLam)
    plot(lambda, -y, type = "b",
      main = "Log marginal density of the ensemble", 
      xlab = "Values of lambda in Box-Cox transformation", ylab = "")
    abline(v = opt$minimum, lty = 2)
    par(op)
  }

  ## invert Box Cox

  "iBoxCox" <- function(y, lambda = 1, tol = 1e-3)
  {
    if (abs(lambda) < tol)
      exp(y)
    else (lambda * y + 1)^(1/lambda)
  }

  ## sample on transformed scale (reuse snippet)

  lambda <- opt$minimum
  EmEB <- BoxCoxSGP(Em0, lambda = lambda)

  sam <- rsamSGP(EmEB, Z = matrix(grid), n = 50)
  sam <- iBoxCox(sam, lambda = lambda)

  if (interactive()) {
    op <- par(ask = TRUE)
    matplot(grid, t(sam), type = "l", 
      main = "Random samples after Box-Cox transformation", xlab = "x",
      ylab = "foo(x)", lty = 1)
    lines(grid, foo(grid), lwd = 2)
    points(drop(Em0$X), Em0$f, pch = 16)
    par(op)
  }

})

#### use mixture of transformations, run with 'eval(SGPexampleMix)'

"SGPexampleMix" <- expression({

  ## check we have run SGPexampleEB

  if (!exists("BoxCox", envir = .GlobalEnv, mode = "function"))
    stop("Must do 'eval(SGPexampleEB)' first")

  ## set up a prior for lambda

  lambda <- seq(from = -3, to = 3, by = 0.5)
  plam <- rep(1, length(lambda)) # doesn't need to be normalised

  ## create a list of emulators, with BoxCox attribute

  if (!exists("myEm", inherits = FALSE))
    stop("Do 'eval(SGPexample)' first")

  Em0 <- myEm
  emList <- lapply(lambda, function(lam) BoxCoxSGP(Em0, lambda = lam))

  ## updated probabilities

  plamNew <- plam * sapply(emList, function(em) exp(-attr(em, "BoxCox")$nls))
  plamNew <- plamNew / sum(plamNew)

  ## summary by sampling

  N <- 50
  ni <- floor(N * plamNew + 0.5) # or use round()

  sam <- matrix(NA, 0, length(grid))
  for (i in which(ni > 0)) {
    em <- emList[[i]]
    lam <- attr(em, "BoxCox")$lambda
    sami <- rsamSGP(em, Z = matrix(grid), n = ni[i])
    sam <- rbind(sam, iBoxCox(sami, lambda = lam))
  }

  ## colour plot by lambda

  if (interactive()) {
    op <- par(ask = TRUE)
    matplot(grid, t(sam), type = "l", 
      main = "Random samples from Box-Cox mixture", xlab = "x",
      ylab = "foo(x)", col = 1 + rep(1:sum(ni > 0), ni[ni > 0]), lty = 1)
    lines(grid, foo(grid), lwd = 2)
    points(drop(Em0$X), Em0$f, pch = 16)

    nonZeros <- paste("p(lam = ", lambda, ") = ",
      round(100 * plamNew), "%", sep = "")[ni > 0]
    legend(x = "bottomright", legend = nonZeros,
      col = 1 + 1:sum(ni > 0), lty = 1, bty = "n")

    par(op)
  }

})