BACCO calibrator - modified for constant mean models and to allow for the 
possiblity of a single physical variable or a single tuning parameter.

## For details, one should refer to the BACCO paper by Hankin, 
## the Kennedy and O'Hagan papers referenced in the BACCO paper
## (especially the technical report) and the BACCO program.


## To begin, make sure you have loaded adapt, mvtnorm, emulator, and 
## calibrator (adapt and mvtnorm should be loaded first)
## using the package manager.

## This document gives some details about setting up the initial
## variables, parameters, and priors.  
## It also briefly describes and provides modified versions of
## several of the functions used in the calibrator package.  The
## modifications are of two general types.  One type is to restrict 
## to the case of the constant mean model.  This mainly involves 
## modifications to the H1 and H2 functions.  The other type is 
## to make the functions work more generally.  Several of the functions 
## are designed for the toy example only.  To generalize them 
## I have added two variables -- Nphys = the number of explanatory
## variables used in the physical model and Ntune = the number of 
## tuning (calibration) variables used in the computer code. 
## These have been added as arguments to many of the functions. 
## I also discovered that when there is only one explanatory variable 
## in the physical model or one tuning variable, many functions 
## return scalars or vectors rather than matrices.  Most of the 
## code seems to assume that objects are matrices.  Thus, in those 
## cases where Nphys = 1 or Ntune = 1, some additional lines of 
## programming (mainly if - else statements) have been added to 
## turn scalars into matrices or modify vectors.


## The next step in running calibrator is to load these modified 
## functions into R.  This can be done by pasting them into 
## R from this document.  Below are the functions that should 
## be pasted in.  At the end of this document I describe a little 
## more about these functions.


phi.fun.ex <- function (rho, lambda, psi1, psi1.apriori, psi2, psi2.apriori, theta.apriori, power, Nphys, Ntune) 
{
    pdm.maker.psi1 <- function(psi1, Nphys, Ntune) {

        if (Nphys == 1) {
            jj.omega_x <- diag(psi1[1], nrow=length(psi1[1]))
            rownames(jj.omega_x) <- names(psi1[1])
            colnames(jj.omega_x) <- names(psi1[1])            
                         }
        else {
            jj.omega_x <- diag(psi1[1:Nphys])
            rownames(jj.omega_x) <- names(psi1[1:Nphys])
            colnames(jj.omega_x) <- names(psi1[1:Nphys])
             }

        if (Ntune == 1) {
             Nt1 <- Nphys + 1
             Nt3 <- Nphys + 2
             jj.omega_t <- diag(psi1[Nt1], nrow=length(psi1[Nt1]))
             rownames(jj.omega_t) <- names(psi1[Nt1])
             colnames(jj.omega_t) <- names(psi1[Nt1])
             sigma1squared <- psi1[Nt3]
                        }

         else {
             Nt1 <- Nphys + 1
             Nt2 <- Nphys + Ntune
             Nt3 <- Nphys + Ntune + 1
             jj.omega_t <- diag(psi1[Nt1:Nt2])
             rownames(jj.omega_t) <- names(psi1[Nt1:Nt2])
             colnames(jj.omega_t) <- names(psi1[Nt1:Nt2])
             sigma1squared <- psi1[Nt3]
               }
        
        return(list(omega_x = jj.omega_x, omega_t = jj.omega_t, 
            sigma1squared = sigma1squared))
    }

    pdm.maker.psi2 <- function(psi2, Nphys, Ntune) {

         if (Nphys == 1) {
             jj.omegastar_x <- diag(psi2[1], nrow=length(psi2[1]))
             sigma2squared <- psi2[2]
 
                         }
         else {
              Np1 <- Nphys + 1
              jj.omegastar_x <- diag(psi2[1:Nphys])
              sigma2squared <- psi2[Np1]
              }

            return(list(omegastar_x = jj.omegastar_x, sigma2squared = 
                sigma2squared))
    }

    jj.mean <- theta.apriori$mean
    jj.V_theta <- theta.apriori$sigma

 if (Ntune == 1) {
    dim(jj.V_theta) = c(1,1)
                 }

    jj.discard.psi1 <- pdm.maker.psi1(psi1, Nphys, Ntune)
    jj.omega_t <- jj.discard.psi1$omega_t
    jj.omega_x <- jj.discard.psi1$omega_x
    jj.sigma1squared <- jj.discard.psi1$sigma1squared
    jj.discard.psi2 <- pdm.maker.psi2(psi2, Nphys, Ntune)
    jj.omegastar_x <- jj.discard.psi2$omegastar_x
    jj.sigma2squared <- jj.discard.psi2$sigma2squared
    jj.omega_t.upper <- chol(jj.omega_t)
    jj.omega_t.lower <- t(jj.omega_t.upper)
    jj.omega_x.upper <- chol(jj.omega_x)
    jj.omega_x.lower <- t(jj.omega_x.upper)
    jj.a <- solve(solve(jj.V_theta) + 2 * jj.omega_t, solve(jj.V_theta, 
        jj.mean))
    jj.b <- t(2 * solve(solve(jj.V_theta) + 2 * jj.omega_t) %*% 
        jj.omega_t)
    jj.c <- jj.sigma1squared/sqrt(det(diag(nrow = nrow(jj.V_theta)) + 
        2 * jj.V_theta %*% jj.omega_t))
    names(jj.c) <- "ht.fun.precalc"
    jj.A <- solve(jj.V_theta + solve(jj.omega_t)/4)
    jj.A.upper <- chol(jj.A)
    jj.A.lower <- t(jj.A.upper)
    list(rho = rho, lambda = lambda, psi1 = psi1, psi1.apriori = psi1.apriori,  psi2 = psi2, psi2.apriori = psi2.apriori, theta.apriori = theta.apriori,  power = power, omega_x = jj.omega_x, omega_t = jj.omega_t, 
omegastar_x = jj.omegastar_x, sigma1squared = jj.sigma1squared, sigma2squared = jj.sigma2squared, omega_x.upper = jj.omega_x.upper,  omega_x.lower = jj.omega_x.lower, omega_t.upper = jj.omega_t.upper,  omega_t.lower = jj.omega_t.lower, a = jj.a, b = jj.b,  c = jj.c, A = jj.A, A.upper = jj.A.upper, A.lower = jj.A.lower)
}



phi.change <- function (phi.fun, old.phi = NULL, rho = NULL, lambda = NULL, 
    psi1 = NULL, psi1.apriori = NULL, psi1.apriori.mean = NULL, 
    psi1.apriori.sigma = NULL, psi2 = NULL, psi2.apriori = NULL, 
    psi2.apriori.mean = NULL, psi2.apriori.sigma = NULL, theta.apriori = NULL, theta.apriori.mean = NULL, theta.apriori.sigma = NULL, power = NULL, 
Nphys, Ntune) 
{
    if (is.null(old.phi)) {
        old.phi <- phi.ex
    }
    if (is.null(rho)) {
        rho <- old.phi$rho
    }
    if (is.null(lambda)) {
        lambda <- old.phi$lambda
    }
    if (is.null(psi1)) {
        psi1 <- old.phi$psi1
    }
    else {
        if (is.null(names(psi1))) {
            names(psi1) <- names(old.phi$psi1)
        }
    }
    if (is.null(psi1.apriori)) {
        psi1.apriori <- old.phi$psi1.apriori
    }
    if (!is.null(psi1.apriori.mean)) {
        psi1.apriori$mean <- psi1.apriori.mean
        if (is.null(names(psi1.apriori.mean))) {
            names(psi1.apriori$mean) <- names(old.phi$psi1.apriori$mean)
        }
    }
    if (!is.null(psi1.apriori.sigma)) {
        psi1.apriori$sigma <- psi1.apriori.sigma
        if (is.null(rownames(psi1.apriori.sigma))) {
            rownames(psi1.apriori$sigma) <- rownames(old.phi$psi1.apriori$sigma)
        }
        if (is.null(colnames(psi1.apriori.sigma))) {
            colnames(psi1.apriori$sigma) <- colnames(old.phi$psi1.apriori$sigma)
        }
    }
    if (is.null(psi2)) {
        psi2 <- old.phi$psi2
    }
    else {
        if (is.null(names(psi2))) {
            names(psi2) <- names(old.phi$psi2)
        }
    }
    if (is.null(psi2.apriori)) {
        psi2.apriori <- old.phi$psi2.apriori
        if (!is.null(psi2.apriori.mean)) {
            psi2.apriori$mean <- psi2.apriori.mean
            if (is.null(names(psi2.apriori.mean))) {
                names(psi2.apriori$mean) <- names(old.phi$psi2.apriori$mean)
            }
        }
        if (!is.null(psi2.apriori.sigma)) {
            psi2.apriori$sigma <- psi2.apriori.sigma
            if (is.null(rownames(psi2.apriori.sigma))) {
                rownames(psi2.apriori$sigma) <- rownames(old.phi$psi2.apriori$sigma)
            }
            if (is.null(colnames(psi2.apriori.sigma))) {
                colnames(psi2.apriori$sigma) <- colnames(old.phi$psi2.apriori$sigma)
            }
        }
    }
    if (is.null(theta.apriori)) {
        theta.apriori <- old.phi$theta.apriori
        if (!is.null(theta.apriori.mean)) {
            theta.apriori$mean <- theta.apriori.mean
            if (is.null(names(theta.apriori.mean))) {
                names(theta.apriori$mean) <- names(old.phi$theta.apriori$mean)
            }
        }
        if (!is.null(theta.apriori.sigma)) {
            theta.apriori$sigma <- theta.apriori.sigma
            if (is.null(rownames(theta.apriori.sigma))) {
                rownames(theta.apriori$sigma) <- rownames(old.phi$theta.apriori$sigma)
            }
            if (is.null(colnames(theta.apriori.sigma))) {
                colnames(theta.apriori$sigma) <- colnames(old.phi$theta.apriori$sigma)
            }
        }
    }
    if (is.null(power)) {
        power <- old.phi$power
    }
    return(phi.fun(rho = rho, lambda = lambda, psi1 = psi1, psi1.apriori = psi1.apriori, psi2 = psi2, psi2.apriori = psi2.apriori, theta.apriori = theta.apriori, power = power, Nphys = Nphys, Ntune = Ntune))
}

H1.ex<-function(D1)
{
out<-rep(1,dim(D1)[1])
dim(out)<-c(dim(D1)[1],1)
return(out)
}


stage1 <- function (D1, y, H1, maxit, method = "Nelder-Mead", directory = ".",  do.filewrite = FALSE, do.print = TRUE, phi.fun, lognormally.distributed = FALSE, include.prior = TRUE, phi, Nphys, Ntune) 
{
    if (do.filewrite & !isTRUE(file.info(directory)$isdir)) {
        stop("do.filewrite = TRUE; directory name supplied does not exist")
    }
    f <- function(candidate) {
        phi.temp <- phi.change(phi.fun = phi.fun, old.phi = phi, 
            psi1 = exp(candidate), Nphys = Nphys, Ntune = Ntune)
        f.out <- p.eqn4.supp(D1 = D1, y = y, H1 = H1, lognormally.distributed = lognormally.distributed, 
            include.prior = include.prior, return.log = TRUE, 
            phi = phi.temp)
        if (do.print) {
            print(candidate)
            print(f.out)
        }
        if (do.filewrite) {
            save(phi.temp, f.out, ascii = TRUE, file = file.path(directory, 
                paste("stage1.", gsub("[ :]", "_", date()), sep = "")))
            save(phi.temp, f.out, ascii = TRUE, file = file.path(directory, 
                "stage1.latest"))
        }
        return(f.out)
    }
    psi1.start <- log(phi$psi1)
    e <- optim(psi1.start, f, method = method, control = list(fnscale = -1, 
        maxit = maxit, trace = 100))
    phi.out <- phi.change(phi.fun = phi.fun, old.phi = phi, psi1 = exp(e$par), Nphys = Nphys, Ntune = Ntune)
    return(invisible(phi.out))
}



H2.ex<-function(D2)
{
out<-rep(1,dim(D2)[1])
dim(out)<-c(dim(D2)[1],1)
return(out)
}


E.theta.ex <- function (D2 = NULL, H1 = NULL, x1 = NULL, x2 = NULL, phi, give.mean = TRUE) 
{
    if (give.mean) {
        m_theta <- phi$theta.apriori$mean
        return(H1(D1.fun(D2, t.vec = m_theta)))
    }
    else {
        out <- 0
        dim(out) = c(1,1)
        return(out)
    }
}

Edash.theta.ex <- function (x, t.vec, k, H1, fast.but.opaque = FALSE, a = NULL, b = NULL, phi = NULL) 
{
    if (fast.but.opaque) {
        edash.mean <- a + crossprod(b, t.vec[k, ])
    }
    else {
        V_theta <- phi$theta.apriori$sigma
        m_theta <- phi$theta.apriori$mean
        omega_t <- phi$omega_t
        edash.mean <- solve(solve(V_theta) + 2 * omega_t, solve(V_theta, 
            m_theta) + 2 * crossprod(omega_t, t.vec[k, ]))
    }
    jj <- as.vector(edash.mean)
    names(jj) <- rownames(edash.mean)
    edash.mean <- jj
    return(H1(D1.fun(x, edash.mean)))
}




extractor.ex <- function(D1, Nphys, Ntune)
{
  Np1 = Nphys + 1
  Npt = Nphys + Ntune
  return(list(x.star=D1[, 1:Nphys, drop=FALSE],
              t.vec=D1[, Np1:Npt, drop=FALSE]))
}

tt.fun <- function (D1, extractor, x.i, x.j, test.for.symmetry = FALSE, 
    method = 1, phi, Nphys, Ntune) 
{
    V_theta <- phi$theta.apriori$sigma
 
if (Ntune == 1) {
    dim(V_theta) <- c(1,1)
                }

    m_theta <- phi$theta.apriori$mean
    x.star <- extractor(D1, Nphys, Ntune)$x.star
    t.vec <- extractor(D1, Nphys, Ntune)$t.vec
    omega_x <- phi$omega_x
    omega_t <- phi$omega_t
    sigma1squared <- phi$sigma1squared
    A.lower <- phi$A.lower
    if (method == 0) {
        x.jframe <- kronecker(rep(1, nrow(t.vec)), t(x.j))
        x.iframe <- kronecker(rep(1, nrow(t.vec)), t(x.i))
        attributes(x.jframe) <- attributes(x.star)
        attributes(x.iframe) <- attributes(x.star)
        exp.x.jk <- exp(-dists.2frames(x.jframe, x.star, omega_x))
        exp.x.il <- exp(-dists.2frames(x.star, x.iframe, omega_x))
        exp.t.kl <- exp(-dists.2frames(t.vec, t.vec, omega_t)/2)
    }
    else {
        exp.x.jk <- kronecker(exp(-dists.2frames(t(x.j), x.star, 
            omega_x)), t(rep(1, nrow(D1))))
        exp.x.il <- kronecker(exp(-t(dists.2frames(x.star, t(x.i), 
            omega_x))), rep(1, nrow(D1)))
        exp.t.kl <- exp(-dists.2frames(t.vec, t.vec, omega_t)/2)
        rownames(exp.x.jk) <- rownames(D1)
        colnames(exp.x.jk) <- rownames(D1)
        rownames(exp.x.il) <- rownames(D1)
        colnames(exp.x.il) <- rownames(D1)
    }
    mtheta.minus.tmean <- matrix(0, nrow(x.star), nrow(x.star))
    for (k in 1:nrow(x.star)) {
        if (test.for.symmetry == TRUE) {
            l.vector <- 1:nrow(x.star)
        }
        else {
            l.vector <- k:nrow(x.star)
        }
        for (l in l.vector) {
            jj <- m_theta - (t.vec[k, ] + t.vec[l, ])/2
            A <- phi$A
            mtheta.minus.tmean[k, l] <- exp(-as.vector(crossprod(jj, 
                A) %*% jj)/2)
        }
    }
    if (test.for.symmetry == FALSE) {
        mtheta.minus.tmean <- symmetrize(mtheta.minus.tmean)
    }
    fish <- 1/sqrt(det(diag(nrow = nrow(V_theta)) + 4 * V_theta %*% 
        omega_t))
    out <- sigma1squared^2 * fish * exp.x.jk * exp.x.il * exp.t.kl * 
        mtheta.minus.tmean
    colnames(out) <- rownames(out)
    return(out)
}

ht.fun <- function (x.i, x.j, D1, extractor, Edash.theta, H1, fast.but.opaque = TRUE,  x.star = NULL, t.vec = NULL, phi, Nphys, Ntune) 
{
    V_theta <- phi$theta.apriori$sigma

if (Ntune == 1) {
    dim(V_theta) <- c(1,1)
                }

    m_theta <- phi$theta.apriori$mean
    omega_x <- phi$omega_x
    omega_t <- phi$omega_t
    if (fast.but.opaque == TRUE) {
        if (is.null(x.star) | is.null(t.vec)) {
            stop("argument 'fast.but.opaque' being TRUE means caller has to supply x.star and t.vec (use extractor.ex() for this)")
        }
    }
    else {
        x.star <- extractor(D1, Nphys, Ntune)$x.star
        t.vec <- extractor(D1, Nphys, Ntune)$t.vec
    }
    f1 <- function(k) {
        as.vector(quad.form(omega_x, x.i - x.star[k, ]))
    }
    matrix.in.f2 <- 2 * V_theta + solve(omega_t)
    f2 <- function(k) {
        as.vector(quad.form.inv(matrix.in.f2, m_theta - t.vec[k, 
            ]))
    }
    jj.a <- phi$a
    jj.b <- phi$b
    jj.c <- phi$c
    f3 <- function(k) {
        if (fast.but.opaque == TRUE) {
            return(Edash.theta(x = x.j, t.vec = t.vec, k = k, 
                H1 = H1, 
                phi = phi))
        }
        else {
            return(Edash.theta(x = x.j, t.vec = t.vec, k = k, 
                H1 = H1, phi = phi))
        }
    }
    out <- sapply(1:nrow(t.vec), function(i) {
        jj.c * exp(-f1(i) - f2(i)) * f3(i)
    })
    dim(out) = c(1,length(out))
    rownames(out) <- colnames(f3(1))
    colnames(out) <- rownames(t.vec)
return(t(out))
}



V.fun <- function (D1, D2, H1, H2, extractor, E.theta, Edash.theta, give.answers = FALSE, test.for.symmetry = FALSE, phi, Nphys, Ntune) 
{
    lambda <- phi$lambda
    rho <- phi$rho
    x.star <- extractor(D1, Nphys, Ntune)$x.star
    t.vec <- extractor(D1, Nphys, Ntune)$t.vec
    v1.d1 <- V1(D1 = D1, other = NULL, phi = phi)
    v1.d1.inv <- solve(v1.d1)
    w1.d1 <- W1(D1 = D1, H1 = H1, phi = phi)
    h1.d1 <- H1(D1)
    w1.blah <- crossprod(w1.d1, crossprod(h1.d1, v1.d1.inv))
    v1.blah <- crossprod(v1.d1.inv, h1.d1) %*% w1.d1
    v1.blah.di.blah <- v1.blah %*% crossprod(h1.d1, v1.d1.inv)
    line1 <- line2 <- line3 <- line4 <- line5 <- line6 <- matrix(0, 
        nrow(D2), nrow(D2))
    tvec.arbitrary <- phi$theta.apriori$mean
    for (i in 1:nrow(D2)) {
        if (test.for.symmetry == FALSE) {
            j.vector <- i:nrow(D2)
        }
        else {
            j.vector <- 1:nrow(D2)
        }
        for (j in j.vector) {
            line1[i, j] <- V1(D1 = D1.fun(x.star = D2[i, ], t.vec = tvec.arbitrary),  other = D1.fun(x.star = D2[j, ], t.vec = tvec.arbitrary), 
phi = phi)
            jj.tt <- tt.fun(D1 = D1, extractor = extractor, x.i = D2[i, 
                ], x.j = D2[j, ], phi = phi, Nphys = Nphys, Ntune = Ntune)
            line2[i, j] <- tr(crossprod(v1.d1.inv, jj.tt))
            jj.hh <- hh.fun(x.i = D2[i, ], x.j = D2[j, ], E.theta = E.theta, 
                H1 = H1, phi = phi)
            line3[i, j] <- tr(crossprod(w1.d1, jj.hh))
            jj.th <- ht.fun(x.i = D2[j, ], x.j = D2[i, ], D1 = D1, 
                extractor = extractor, Edash.theta = Edash.theta, 
                H1 = H1, fast.but.opaque = TRUE, x.star = x.star, 
                t.vec = t.vec, phi = phi, Nphys = Nphys, Ntune = Ntune)
            line4[i, j] <- tr(w1.blah %*% jj.th)
            jj.ht <- ht.fun(x.i = D2[i, ], x.j = D2[j, ], D1 = D1, 
                extractor = extractor, Edash.theta = Edash.theta, 
                H1 = H1, fast.but.opaque = TRUE, x.star = x.star, 
                t.vec = t.vec, phi = phi, Nphys = Nphys, Ntune = Ntune)
            line5[i, j] <- tr(crossprod(v1.blah, jj.ht))
            line6[i, j] <- tr(crossprod(jj.tt, v1.blah.di.blah))
        }
    }
    C.m <- +line1 - line2 + line3 - line4 - line5 + line6
    if (test.for.symmetry == FALSE) {
        C.m <- symmetrize(C.m)
    }
    C.m <- rho^2 * C.m
    rownames(C.m) <- rownames(D2)
    colnames(C.m) <- rownames(D2)
    C.lambda <- lambda * diag(nrow = nrow(D2))
    C.v2 <- V2(D2, other = NULL, phi = phi)
    out <- C.lambda + C.v2 + C.m

    if (give.answers) {
        attributes(line1) <- attributes(out)
        attributes(line2) <- attributes(out)
        attributes(line3) <- attributes(out)
        attributes(line4) <- attributes(out)
        attributes(line5) <- attributes(out)
        attributes(line6) <- attributes(out)
        return(list(line1 = line1, line2 = line2, line3 = line3, 
            line4 = line4, line5 = line5, line6 = line6, v1.d1 = v1.d1, 
            v1.d1.inv = v1.d1.inv, w1.d1 = w1.d1, h1.d1 = h1.d1, 
            w1.blah = w1.blah, v1.blah = v1.blah, v1.blah.di.blah = v1.blah.di.blah, 
            C.m = C.m, C.lambda = C.lambda, C.v2 = C.v2, out = out))
    }
    else {
        return(out)
    }
}



W2 <- function (D2, H2, V, det = FALSE, Ntune) 
{
    out <- quad.form(solve(V), H2(D2))
if (length(out) == 1) {
                     dim(out) <-c(1,1)
                     }
if (Ntune == 1) {
                 dim(out) <-c(1,1)
                }
    if (det) {
        return(1/det(out))
    }
    else {
        return(solve(out))
    }
}


beta2hat.fun <- function (D1, D2, H1, H2, V = NULL, z, etahat.d2, extractor, E.theta, Edash.theta, phi, Nphys, Ntune) 
{
    rho <- phi$rho
    if (is.null(V)) {
        V <- V.fun(D1 = D1, D2 = D2, H1 = H1, H2 = H2, extractor = extractor, E.theta = E.theta, Edash.theta = Edash.theta, phi = phi, Nphys, Ntune = Ntune)
    }
if (Ntune == 1) {
                 out <- crossprod(W2(D2, H2, V, Ntune = Ntune), 
                 crossprod(H2(D2),
                 solve(V, t(z - rho * etahat.d2))))
                 return(out)
                }
out <- crossprod(W2(D2, H2, V, Ntune = Ntune), crossprod(H2(D2), 
          solve(V, z - rho * etahat.d2)))
return(out)
}


t.fun <- function (x, D1, extractor, phi, Nphys, Ntune) 
{
    sigma1squared <- phi$sigma1squared
    rho <- phi$rho
    omega_x <- phi$omega_x
    omega_t <- phi$omega_t
    x.star <- extractor(D1, Nphys, Ntune)$x.star
    t.vec <- extractor(D1, Nphys, Ntune)$t.vec
    m_theta <- phi$theta.apriori$mean
    V_theta <- phi$theta.apriori$sigma
    jj.mat <- solve(2 * V_theta + solve(omega_t))
    prod.1 <- sigma1squared/sqrt(det(diag(nrow = nrow(omega_t)) + 
        2 * crossprod(V_theta, omega_t)))
    out <- rep(NA, nrow(D1))
    for (j in 1:length(out)) {
        jj.x <- as.vector(x.star[j, ] - x)
        jj.m <- as.vector(m_theta - t.vec[j, ])
        prod.2 <- exp(-quad.form(omega_x, jj.x))
        prod.3 <- exp(-quad.form(jj.mat, jj.m))
        out[j] <- prod.1 * prod.2 * prod.3
    }
    names(out) <- rownames(t.vec)
    return(out)
}



etahat <- function (D1, D2, H1, y, E.theta, extractor, phi, Nphys, Ntune) 
{
    jj.hfun <- E.theta(D2 = D2, H1 = H1, phi = phi, give.mean = TRUE)
    jj.beta <- beta1hat.fun(D1 = D1, H1 = H1, y = y, phi = phi)
    bit1 <- drop(jj.hfun %*% jj.beta)
    b <- beta1hat.fun(D1 = D1, H1 = H1, y = y, phi = phi)
    f <- function(x) {
        t.fun(x, D1 = D1, extractor = extractor, phi = phi, Nphys = Nphys, Ntune = Ntune)
    }
    jj.tfun <- apply(D2, 1, f)
    bit2 <- crossprod(jj.tfun, solve(V1(D1, phi = phi), y - H1(D1) %*% 
        b))
    return(drop(bit1 + bit2))
}



p.page4 <- function (D1, D2, H1, H2, V, y, z, E.theta, Edash.theta, extractor, include.prior = FALSE, lognormally.distributed = FALSE, return.log = FALSE, phi, Nphys, Ntune) 
{
    if (is.null(V)) {
        V <- V.fun(D1 = D1, D2 = D2, H1 = H1, H2 = H2, extractor = extractor, E.theta = E.theta, Edash.theta = Edash.theta, give.answers = FALSE, phi = phi, Nphys = Nphys, Ntune = Ntune)
    }
    etahat.d2 <- etahat(D1 = D1, D2 = D2, H1 = H1, y = y, E.theta = E.theta, extractor = extractor, phi = phi, Nphys, Ntune)
    beta2hat <- beta2hat.fun(D1 = D1, D2 = D2, H1 = H1, H2 = H2, 
        V = V, z = z, etahat.d2 = etahat.d2, extractor = extractor, 
        E.theta = E.theta, Edash.theta = Edash.theta, phi = phi, Nphys = Nphys, Ntune = Ntune)
    h2.d2 <- H2(D2)
    rho <- phi$rho
if (Ntune == 1) {
                 vec <-t(z) - h2.d2 %*% beta2hat - rho * etahat.d2
                }
else            {
                 vec <- z - h2.d2 %*% beta2hat - rho * etahat.d2
                } 

    if (return.log) {
        bit1 <- log(W2(D2, H2, V, det = TRUE, Ntune = Ntune))/2 - 
                       log(det(V))/2
        bit2 <- -0.5 * quad.form.inv(V, vec)
        if (include.prior) {
            return(bit1 + bit2 + log(prob.psi2(phi, lognormally.distributed = lognormally.distributed)))
        }
        else {
            return(bit1 + bit2)
        }
    }
    else {
        bit1 <- sqrt(W2(D2, H2, V, det = TRUE, Ntune = Ntune)/det(V))
        bit2 <- exp(-0.5 * quad.form.inv(V, vec))
        if (include.prior) {
            return(bit1 * bit2 * prob.psi2(phi, lognormally.distributed = lognormally.distributed))
        }
        else {
            return(bit1 * bit2)
        }
    }
}



stage2 <- function (D1, D2, H1, H2, y, z, maxit, method = "Nelder-Mead", 
    directory = ".", do.filewrite = FALSE, do.print = TRUE, extractor, 
    phi.fun, E.theta, Edash.theta, lognormally.distributed = FALSE, 
    use.standin = FALSE, phi, Nphys, Ntune) 
{
    if (do.filewrite & !isTRUE(file.info(directory)$isdir)) {
        stop("do.filewrite = TRUE; directory name supplied does not exist")
    }

Np3 = Nphys + 3

    f <- function(candidate) {
        phi.temp <- phi.change(phi.fun = phi.fun, old.phi = phi, 
            rho = exp(candidate[1]), lambda = exp(candidate[2]), 
            psi2 = exp(candidate[3:Np3]), Nphys = Nphys, Ntune = 
            Ntune)
        if (use.standin) {
            V.temp <- 1 + diag(3, nrow = nrow(D2))
        }
        else {
            V.temp <- V.fun(D1 = D1, D2 = D2, H1 = H1, H2 = H2, 
                extractor = extractor, E.theta = E.theta, Edash.theta = Edash.theta,  phi = phi.temp, Nphys = Nphys, Ntune = Ntune)
        }
        f.out <- drop(p.page4(D1 = D1, D2 = D2, H1 = H1, H2 = H2, 
            V = V.temp, z = z, y = y, E.theta = E.theta, Edash.theta = Edash.theta,  extractor = extractor, include.prior = FALSE, lognormally.distributed = lognormally.distributed,  return.log = TRUE, 
phi = phi.temp, Nphys = Nphys, Ntune = Ntune))
        if (do.print) {
            print(candidate)
            print(f.out)
        }
        if (do.filewrite) {
            save(phi.temp, f.out, ascii = TRUE, file = file.path(directory, 
                paste("stage2.", gsub("[ :]", "_", date()), sep = "")))
            save(phi.temp, f.out, ascii = TRUE, file = file.path(directory, 
                "stage2.latest"))
        }
        return(f.out)
    }
    rho.lambda.psi2.start <- log(c(phi$rho, phi$lambda, phi$psi2))
    e <- optim(rho.lambda.psi2.start, f, method = method, control = list(fnscale = -1, 
        trace = 100, maxit = maxit))
    jj <- exp(e$par)
    phi.out <- phi.change(phi.fun = phi.fun, old.phi = phi, rho = jj[1], 
        lambda = jj[2], psi2 = jj[3:Np3], Nphys = Nphys, Ntune = Ntune)
    return(invisible(phi.out))
}



stage3 <- function (D1, D2, H1, H2, d, maxit, trace = 100, method = "Nelder-Mead", 
    directory = ".", do.filewrite = FALSE, do.print = TRUE, include.prior = TRUE, 
    lognormally.distributed = FALSE, theta.start = NULL, phi) 
{

p.eqn8.supp <- function (theta, D1, D2, H1, H2, d, include.prior = FALSE, lognormally.distributed = FALSE, 
    return.log = FALSE, phi) 
{

p.eqn8.supp.vector <- function (theta, D1, D2, H1, H2, d, include.prior = FALSE, lognormally.distributed = FALSE, 
    return.log = FALSE, phi) 
{
p.theta <- function (theta, phi, lognormally.distributed = FALSE) 
{

if (length(theta) == 1) {
             return(dnorm(x = theta, mean = phi$theta.aprior$mean,
                    sd = sqrt(phi$theta.aprior$sigma)))
                        }
else   {
    if (lognormally.distributed) {
        if (any(x) < 0) {
            return(0)
                        }
        return(dmvnorm(x = log(theta), mean = log(phi$theta.aprior$mean), 
            sigma = phi$theta.aprior$sigma))
                                 }
       else {
        return(dmvnorm(x = theta, mean = phi$theta.aprior$mean, 
            sigma = phi$theta.aprior$sigma))
             }
          }
}


    betahat <- betahat.fun.koh(theta = theta, D1 = D1, D2 = D2, 
        H1 = H1, H2 = H2, phi = phi, d = d)
    Vd.theta <- Vd(theta = theta, D1 = D1, D2 = D2, phi = phi)
    det.W.theta <- W(theta = theta, D1 = D1, D2 = D2, H1 = H1, 
        H2 = H2, det = TRUE, phi = phi)
    H.theta <- H.fun(theta = theta, D1 = D1, D2 = D2, H1 = H1, 
        H2 = H2, phi = phi)
    mismatch <- d - H.theta %*% betahat
    if (return.log) {
        out <- -log(det(Vd.theta))/2 + log(det.W.theta)/2 - quad.form.inv(Vd.theta, 
            mismatch)/2
        if (include.prior) {
            out <- out + log(p.theta(theta, phi, lognormally.distributed = lognormally.distributed))
        }
    }
    else {
        out <- det(Vd.theta)^(-1/2) * det.W.theta^(1/2) * exp(-0.5 * 
            quad.form.inv(Vd.theta, mismatch))
        if (include.prior) {
            out <- out * p.theta(theta, phi, lognormally.distributed = lognormally.distributed)
        }
    }
    return(as.vector(out))
}


    if (is.vector(theta)) {
        return(p.eqn8.supp.vector(theta = theta, D1 = D1, D2 = D2, 
            H1 = H1, H2 = H2, d = d, include.prior = include.prior, 
            lognormally.distributed = lognormally.distributed, 
            return.log = return.log, phi = phi))
    }

else{
   return(0)
    }
}


    if (do.filewrite & !isTRUE(file.info(directory)$isdir)) {
        stop("do.filewrite = TRUE; directory name supplied does not exist")
    }
    f <- function(theta) {
        print(theta)
        f.out <- p.eqn8.supp(theta, D1 = D1, D2 = D2, H1 = H1, 
            H2 = H2, d = d, include.prior = include.prior, lognormally.distributed = lognormally.distributed, 
            return.log = TRUE, phi = phi)
        if (do.print) {
            print(theta)
            print(f.out)
        }
        if (do.filewrite) {
            save(theta, f.out, ascii = TRUE, file = file.path(directory, 
                paste("stage3.", gsub("[ :]", "_", date()), sep = "")))
            save(theta, f.out, ascii = TRUE, file = file.path(directory, 
                "stage3.latest"))
        }
        return(f.out)
    }
    if (is.null(theta.start)) {
        theta.start <- phi$theta.apriori$mean
    }
    e <- optim(theta.start, f, method = method, control = list(fnscale = -1, 
        maxit = maxit, trace = trace))
    optimal.theta <- e$par
    names(optimal.theta) <- names(phi$theta.apriori$mean)
    return(optimal.theta)
}



## After pasting the functions into R, you must next provide the 
## variables or training data.
## In BACCO calibrator, there are two sets of explanatory variables 
## (training data) that must be defined.  Both sets must be in a matrix.
## Each row of these matrices corresponds to a different observation.

## One set of  training data consists of the values of the explanatory 
## variables for the computer code.
## These explanatory variables include both the actual physical 
## variables and the tuning variables.
## The following are for a one-dimensional example with one 
## physical variable and one tuning variable.  At the end of 
## this document all these definitions are grouped together 
## so they can be easily cut and pasted into R if you want to try them
## out.  I also give a second example (essentially the toy example
## in BACCO) that can also be cut and pasted into R to try the code
## out and a multidimensional example with two physical and three 
## tuning variables.

xcode.training <- latin.hypercube(11, 2)

## The other set of training data consists of the values of the explanatory 
## variables for the physical experiment.  
## These explanatory variables include just the actual physical variables.
## These need to be in a column vector.

xphys.training <- seq(0,1,0.1) 
dim(xphys.training) <- c(11,1)

## It is useful to specify the number of physical variables 
## and the number of tuning parameters.
## In this example, both are 1.

Numphys <- 1
Numtuning <- 1

## For this simple example, we generate the computer code data 
## according to a known function.

codedatagenerator <-function(x, t)
            exp(-1.4*t*x)*cos((3.5*pi*x)) 

## and we assume the physical experiment produces data according 
## to the function (with random error)

physdatagenerator <-function(x)
            exp(-1.4*x)*cos((3.5*pi*x)) + rnorm(1, mean=0, sd=0.1)

## This function assumes the true value of the calibration 
## parameter is 1.

## In BACCO, there are three responses (training values) that must 
## be specified.  All must be row vectors (or matrices with 1 row).
## The first response is the output from the code.

ycode.training <- codedatagenerator(xcode.training[ ,1], 
    xcode.training[ ,2])


## The second response is the observations from the physical experiment.

yphys.training <- physdatagenerator(xphys.training)
dim(yphys.training) <- c(1,11)

## The third response is the vector of all the training data responses,
## which consists of the ycode.training 
## values followed by the yphys.training values.

y.training <- c(ycode.training, yphys.training)

## In order to proceed, one must specify the values of hyperparameters
## and priors used in the  model.  See the paper by Hankin for meanings.


## First are the (starting?) values of the scales for the correlation 
## function for the code plus the variance sigma1squared.  
## Recall that the explanatory variables in the 
## code consist of the physical variables and the tuning variables.
## The first variables listed are the scales for the physical variables
## followed by the scales for the tuning variables, followed by the 
## value for sigma1squared.

psi1.ex <- c(2,2,2)

## Next is the prior distribution for psi1.  This is a two element
## list with the first element the prior mean and the second the prior 
## covariance (sigma) matrix.  Note that actually in this prior 
## one is specifying the prior for the scales and for the variance 
## sigma1squared.
## Thus, in specifying the prior mean, the first entries are the 
## prior means for the psi1 parameters and the last is that for 
## sigma1squared.  Likewise, in specifying the prior sigma, the
## entries correspond to those in the prior mean (first correspond
## to psi1 and last or lower right entry is that for sigma1sqaured.
## Note that this prior is for the logarithms of the parameters.  
## The logarithms of the parameters are assumed to be multivariate
## normal with the specified mean and covariance.

psi1.apriori.ex <- list(mean=c(0,0,1), sigma=1+diag(3))

## Next, are two (starting values?) univatiate parameters rho and lambda.

rho.ex <- 1
lambda.ex <- 1


## Next, are the (starting?) values of the scales for the correlation 
## function for the physical experiment plus the variance sigma2sqaured.  
## Recall that the explanatory variables 
## in the physical experiment consist only of the physical variables.
## Thus, the first variables listed are the scales for the physical 
## variables followed by the value for sigma2squared.

psi2.ex <- c(10,2)

## Next is the prior distribution for psi2.  This is a two element
## list with the first element the prior mean and the second the prior 
## covariance (sigma) matrix.  Note that actually in this prior 
## one is specifying the prior for rho, lambda, the scale parameters, 
## and the variance sigma2squared.
## Thus, in specifying the prior mean, the first two
## entries are the means for rho and lamda.  The next are the
## prior means for the scale parameters and the last is that for 
## sigma2squared.  Likewise, in specifying the prior sigma, the entries
## correspond with those in the prior mean (first two diagonal entries
## correspond to rho and lambda, next entries to the scale 
## parameters, and last to sigma2squared.
## Note that this prior is for the logarithms of the parameters.  
## The logarithms of the parameters are assumed to be multivariate
## normal with the specified mean and covariance.

psi2.apriori.ex <- list(mean=c(0,0,0,1), sigma=0.5+diag(4))

## Next is the "true" (starting?) value of the tuning parmaters.
## It seems to be the case that theta can either be entered 
## as a number or vector of form c(.,.,...).

theta.ex <- 1

## Next is the prior for theta. This is a two element
## list with the first element the prior mean and the second the prior 
## covariance (sigma) matrix. 
## Note that, unlike the previous priors, here we assume the prior for
## theta itself is multivariate normal with the specified mean 
## and covariance.

theta.apriori.ex <- list(mean=0, sigma=5)

## Finally we specify the power for the correlation function.  Unless
## there is reason to believe the models are not smooth, use 2.

power.ex <- 2


## Now we are ready to actually do the calibration.  We first create 
## the phi object (here I call it phi.ex).  
## Then we carry out the stage1, stage2, and stage3 optimizations.
## The stage3 optimization produces the estimate of the tuning 
## (calibration) parameters theta.

phi.ex <- phi.fun.ex(rho.ex, lambda.ex, psi1.ex, psi1.apriori.ex, psi2.ex, psi2.apriori.ex, theta.apriori.ex, power.ex, Numphys, Numtuning)

phi.ex.stage1 <- stage1(D1=xcode.training, y=ycode.training, H1=H1.ex,
  maxit=5, phi.fun=phi.fun.ex, phi=phi.ex, Nphys = Numphys, Ntune = Numtuning)


phi.ex.stage2 <- stage2(D1=xcode.training, D2=xphys.training, H1=H1.ex, H2=H2.ex, y=ycode.training, z=yphys.training, extractor=extractor.ex,
phi.fun=phi.fun.ex, E.theta=E.theta.ex, Edash.theta=Edash.theta.ex,
maxit=3, phi=phi.ex.stage1,  Nphys = Numphys, Ntune = Numtuning)

theta.stage3 <- stage3(D1 = xcode.training, D2 = xphys.training, H1=H1.ex, H2=H2.ex, d = y.training, maxit = 10, phi = phi.ex.stage2)

## Finally, you can use the function p.eqn8.supp to compute 
## the posterior probabilities for values of theta given 
## the data and the values of the scale parameters (psi1). 
## For example, the following calculates the posterior 
## prob at theta.ex, our starting value.

p.eqn8.supp(theta=theta.ex, D1=xcode.training, D2=xphys.training, 
H1=H1.ex, H2=H2.ex, d=y.training, include.prior = TRUE, lognormally.distributed = FALSE, phi=phi.ex.stage2)


## See the help for calibrator to learn more about p.eqn8.supp.


TWO EXAMPLES

## In both of the following examples, one first must load 
## adapt, mvtnorm, emulator, and calibrator into R.
## Then one must paste all the modified functions into 
## R as described previously.
## Note that the priors used in these examples 
## were not chosen for any particular reason.

## EXAMPLE 1 - A 1-dimensional example (as above)


## The data, parameter starting values, and priors.

xcode.training <- latin.hypercube(11, 2)
xphys.training <- seq(0,1,0.1) 
dim(xphys.training) <- c(11,1)
Numphys <- 1
Numtuning <- 1
codedatagenerator <-function(x, t)
            exp(-1.4*t*x)*cos((3.5*pi*x)) 
physdatagenerator <-function(x)
            exp(-1.4*x)*cos((3.5*pi*x)) + rnorm(1, mean=0, sd=0.1)
ycode.training <- codedatagenerator(xcode.training[ ,1], 
    xcode.training[ ,2])
yphys.training <- physdatagenerator(xphys.training)
dim(yphys.training) <- c(1,11)
y.training <- c(ycode.training, yphys.training)


psi1.ex <- c(2,2,2)
psi1.apriori.ex <- list(mean=c(0,0,1), sigma=1+diag(3))
rho.ex <- 1
lambda.ex <- 1
psi2.ex <- c(10,2)
psi2.apriori.ex <- list(mean=c(0,0,0,1), sigma=0.5+diag(4))
theta.ex <- 1
theta.apriori.ex <- list(mean=0, sigma=5)
power.ex <- 2


## Create phi and run the optimizations

phi.ex <- phi.fun.ex(rho.ex, lambda.ex, psi1.ex, psi1.apriori.ex, psi2.ex, psi2.apriori.ex, theta.apriori.ex, power.ex, Numphys, Numtuning)

phi.ex.stage1 <- stage1(D1=xcode.training, y=ycode.training, H1=H1.ex,
  maxit=5, phi.fun=phi.fun.ex, phi=phi.ex, Nphys = Numphys, Ntune = 
  Numtuning)


phi.ex.stage2 <- stage2(D1=xcode.training, D2=xphys.training, H1=H1.ex, H2=H2.ex, y=ycode.training, z=yphys.training, extractor=extractor.ex,
phi.fun=phi.fun.ex, E.theta=E.theta.ex, Edash.theta=Edash.theta.ex,
maxit=3, phi=phi.ex.stage1,  Nphys = Numphys, Ntune = Numtuning)

theta.stage3 <- stage3(D1 = xcode.training, D2 = xphys.training, H1=H1.ex, H2=H2.ex, d = y.training, maxit = 10, phi = phi.ex.stage2)

p.eqn8.supp(theta=theta.ex, D1=xcode.training, D2=xphys.training, 
H1=H1.ex, H2=H2.ex, d=y.training, include.prior = TRUE, phi=phi.ex.stage2)


EXAMPLE 2 - The example with the toy data from Hankin paper

## The data, parameter starting values, and priors.

data(toys)

xcode.training <- D1.toy
xphys.training <- D2.toy 
Numphys <- 2
Numtuning <- 3
ycode.training <- y.toy
yphys.training <-z.toy
y.training <- d.toy


psi1.ex <- structure(c(1.1, 1.2, 1.3, 1.4,1.5,0.7), .Names = c("x", "y",
                      "A","B", "C","s1sq"))
psi1.apriori.ex <- list(mean=rep(0,6),sigma=0.4+diag(6))
rho.ex <- 1
lambda.ex <- 1
psi2.ex <- structure(c(2.1, 2.2, 1), .Names = c("x","y","s2sq"))
psi2.apriori.ex <- list(mean=rep(0,5),sigma=0.2+diag(5))
theta.ex <- theta.toy
theta.apriori.ex <- list(mean=0.1+(1:3),sigma=2.1+diag(3))
power.ex <- 2


## Create phi and run the optimizations

phi.ex <- phi.fun.ex(rho.ex, lambda.ex, psi1.ex, psi1.apriori.ex, psi2.ex, psi2.apriori.ex, theta.apriori.ex, power.ex, Numphys, Numtuning)

phi.ex.stage1 <- stage1(D1=xcode.training, y=ycode.training, H1=H1.ex,
  maxit=5, phi.fun=phi.fun.ex, phi=phi.ex, Nphys = Numphys, Ntune = Numtuning)


phi.ex.stage2 <- stage2(D1=xcode.training, D2=xphys.training, H1=H1.ex, H2=H2.ex, y=ycode.training, z=yphys.training, extractor=extractor.ex,
phi.fun=phi.fun.ex, E.theta=E.theta.ex, Edash.theta=Edash.theta.ex,
maxit=3, phi=phi.ex.stage1,  Nphys = Numphys, Ntune = Numtuning)

theta.stage3 <- stage3(D1 = xcode.training, D2 = xphys.training, H1=H1.ex, H2=H2.ex, d = y.training, maxit = 10, phi = phi.ex.stage2)

p.eqn8.supp(theta=theta.ex, D1=xcode.training, D2=xphys.training, 
H1=H1.ex, H2=H2.ex, d=y.training, include.prior = TRUE, phi=phi.ex.stage2)


A BIT MORE ON THE MODIFIED FUNCTIONS

## A function (denoted phi.fun in BACCO) must be defined 
## that generates some additional parameters from the above. 
## These are then assigned to an object "phi".
## The function phi.fun.toy in BACCO doesn't work if 
## there is only a single physical variable or a single 
## tuning parameter.  The following is a modified version 
## that works in all cases.

phi.fun.ex <- function (rho, lambda, psi1, psi1.apriori, psi2, psi2.apriori, theta.apriori, power, Nphys, Ntune) 
{
    pdm.maker.psi1 <- function(psi1, Nphys, Ntune) {

        if (Nphys == 1) {
            jj.omega_x <- diag(psi1[1], nrow=length(psi1[1]))
            rownames(jj.omega_x) <- names(psi1[1])
            colnames(jj.omega_x) <- names(psi1[1])            
                         }
        else {
            jj.omega_x <- diag(psi1[1:Nphys])
            rownames(jj.omega_x) <- names(psi1[1:Nphys])
            colnames(jj.omega_x) <- names(psi1[1:Nphys])
             }

        if (Ntune == 1) {
             Nt1 <- Nphys + 1
             Nt3 <- Nphys + 2
             jj.omega_t <- diag(psi1[Nt1], nrow=length(psi1[Nt1]))
             rownames(jj.omega_t) <- names(psi1[Nt1])
             colnames(jj.omega_t) <- names(psi1[Nt1])
             sigma1squared <- psi1[Nt3]
                        }

         else {
             Nt1 <- Nphys + 1
             Nt2 <- Nphys + Ntune
             Nt3 <- Nphys + Ntune + 1
             jj.omega_t <- diag(psi1[Nt1:Nt2])
             rownames(jj.omega_t) <- names(psi1[Nt1:Nt2])
             colnames(jj.omega_t) <- names(psi1[Nt1:Nt2])
             sigma1squared <- psi1[Nt3]
               }
        
        return(list(omega_x = jj.omega_x, omega_t = jj.omega_t, 
            sigma1squared = sigma1squared))
    }

    pdm.maker.psi2 <- function(psi2, Nphys, Ntune) {

         if (Nphys == 1) {
             jj.omegastar_x <- diag(psi2[1], nrow=length(psi2[1]))
             sigma2squared <- psi2[2]
 
                         }
         else {
              Np1 <- Nphys + 1
              jj.omegastar_x <- diag(psi2[1:Nphys])
              sigma2squared <- psi2[Np1]
              }

            return(list(omegastar_x = jj.omegastar_x, sigma2squared = 
                sigma2squared))
    }

    jj.mean <- theta.apriori$mean
    jj.V_theta <- theta.apriori$sigma

 if (Ntune == 1) {
    dim(jj.V_theta) = c(1,1)
                 }

    jj.discard.psi1 <- pdm.maker.psi1(psi1, Nphys, Ntune)
    jj.omega_t <- jj.discard.psi1$omega_t
    jj.omega_x <- jj.discard.psi1$omega_x
    jj.sigma1squared <- jj.discard.psi1$sigma1squared
    jj.discard.psi2 <- pdm.maker.psi2(psi2, Nphys, Ntune)
    jj.omegastar_x <- jj.discard.psi2$omegastar_x
    jj.sigma2squared <- jj.discard.psi2$sigma2squared
    jj.omega_t.upper <- chol(jj.omega_t)
    jj.omega_t.lower <- t(jj.omega_t.upper)
    jj.omega_x.upper <- chol(jj.omega_x)
    jj.omega_x.lower <- t(jj.omega_x.upper)
    jj.a <- solve(solve(jj.V_theta) + 2 * jj.omega_t, solve(jj.V_theta, 
        jj.mean))
    jj.b <- t(2 * solve(solve(jj.V_theta) + 2 * jj.omega_t) %*% 
        jj.omega_t)
    jj.c <- jj.sigma1squared/sqrt(det(diag(nrow = nrow(jj.V_theta)) + 
        2 * jj.V_theta %*% jj.omega_t))
    names(jj.c) <- "ht.fun.precalc"
    jj.A <- solve(jj.V_theta + solve(jj.omega_t)/4)
    jj.A.upper <- chol(jj.A)
    jj.A.lower <- t(jj.A.upper)
    list(rho = rho, lambda = lambda, psi1 = psi1, psi1.apriori = psi1.apriori,  psi2 = psi2, psi2.apriori = psi2.apriori, theta.apriori = theta.apriori,  power = power, omega_x = jj.omega_x, omega_t = jj.omega_t, 
omegastar_x = jj.omegastar_x, sigma1squared = jj.sigma1squared, sigma2squared = jj.sigma2squared, omega_x.upper = jj.omega_x.upper,  omega_x.lower = jj.omega_x.lower, omega_t.upper = jj.omega_t.upper,  omega_t.lower = jj.omega_t.lower, a = jj.a, b = jj.b,  c = jj.c, A = jj.A, A.upper = jj.A.upper, A.lower = jj.A.lower)
}



## Now we must modify phi.change

phi.change <- function (phi.fun, old.phi = NULL, rho = NULL, lambda = NULL, 
    psi1 = NULL, psi1.apriori = NULL, psi1.apriori.mean = NULL, 
    psi1.apriori.sigma = NULL, psi2 = NULL, psi2.apriori = NULL, 
    psi2.apriori.mean = NULL, psi2.apriori.sigma = NULL, theta.apriori = NULL, theta.apriori.mean = NULL, theta.apriori.sigma = NULL, power = NULL, 
Nphys, Ntune) 
{
    if (is.null(old.phi)) {
        old.phi <- phi.ex
    }
    if (is.null(rho)) {
        rho <- old.phi$rho
    }
    if (is.null(lambda)) {
        lambda <- old.phi$lambda
    }
    if (is.null(psi1)) {
        psi1 <- old.phi$psi1
    }
    else {
        if (is.null(names(psi1))) {
            names(psi1) <- names(old.phi$psi1)
        }
    }
    if (is.null(psi1.apriori)) {
        psi1.apriori <- old.phi$psi1.apriori
    }
    if (!is.null(psi1.apriori.mean)) {
        psi1.apriori$mean <- psi1.apriori.mean
        if (is.null(names(psi1.apriori.mean))) {
            names(psi1.apriori$mean) <- names(old.phi$psi1.apriori$mean)
        }
    }
    if (!is.null(psi1.apriori.sigma)) {
        psi1.apriori$sigma <- psi1.apriori.sigma
        if (is.null(rownames(psi1.apriori.sigma))) {
            rownames(psi1.apriori$sigma) <- rownames(old.phi$psi1.apriori$sigma)
        }
        if (is.null(colnames(psi1.apriori.sigma))) {
            colnames(psi1.apriori$sigma) <- colnames(old.phi$psi1.apriori$sigma)
        }
    }
    if (is.null(psi2)) {
        psi2 <- old.phi$psi2
    }
    else {
        if (is.null(names(psi2))) {
            names(psi2) <- names(old.phi$psi2)
        }
    }
    if (is.null(psi2.apriori)) {
        psi2.apriori <- old.phi$psi2.apriori
        if (!is.null(psi2.apriori.mean)) {
            psi2.apriori$mean <- psi2.apriori.mean
            if (is.null(names(psi2.apriori.mean))) {
                names(psi2.apriori$mean) <- names(old.phi$psi2.apriori$mean)
            }
        }
        if (!is.null(psi2.apriori.sigma)) {
            psi2.apriori$sigma <- psi2.apriori.sigma
            if (is.null(rownames(psi2.apriori.sigma))) {
                rownames(psi2.apriori$sigma) <- rownames(old.phi$psi2.apriori$sigma)
            }
            if (is.null(colnames(psi2.apriori.sigma))) {
                colnames(psi2.apriori$sigma) <- colnames(old.phi$psi2.apriori$sigma)
            }
        }
    }
    if (is.null(theta.apriori)) {
        theta.apriori <- old.phi$theta.apriori
        if (!is.null(theta.apriori.mean)) {
            theta.apriori$mean <- theta.apriori.mean
            if (is.null(names(theta.apriori.mean))) {
                names(theta.apriori$mean) <- names(old.phi$theta.apriori$mean)
            }
        }
        if (!is.null(theta.apriori.sigma)) {
            theta.apriori$sigma <- theta.apriori.sigma
            if (is.null(rownames(theta.apriori.sigma))) {
                rownames(theta.apriori$sigma) <- rownames(old.phi$theta.apriori$sigma)
            }
            if (is.null(colnames(theta.apriori.sigma))) {
                colnames(theta.apriori$sigma) <- colnames(old.phi$theta.apriori$sigma)
            }
        }
    }
    if (is.null(power)) {
        power <- old.phi$power
    }
    return(phi.fun(rho = rho, lambda = lambda, psi1 = psi1, psi1.apriori = psi1.apriori, psi2 = psi2, psi2.apriori = psi2.apriori, theta.apriori = theta.apriori, power = power, Nphys = Nphys, Ntune = Ntune))
}


## To do the stage1 optimization, we need to 
## define the H1 function which returns the column vector of 
## the regression function values.  For standard computer
## experiment models, this will be a column vector of 1s
## and the number of rows in the column will be equal to the 
## number of observations on the code.  
## Note that this is a modification of the H1.toy function in 
## BACCO, which allowed for linear and other regression terms.

H1.ex<-function(D1)
{
out<-rep(1,dim(D1)[1])
dim(out)<-c(dim(D1)[1],1)
return(out)
}

## Now we must modify the function stage1

stage1 <- function (D1, y, H1, maxit, method = "Nelder-Mead", directory = ".",  do.filewrite = FALSE, do.print = TRUE, phi.fun, lognormally.distributed = FALSE, include.prior = TRUE, phi, Nphys, Ntune) 
{
    if (do.filewrite & !isTRUE(file.info(directory)$isdir)) {
        stop("do.filewrite = TRUE; directory name supplied does not exist")
    }
    f <- function(candidate) {
        phi.temp <- phi.change(phi.fun = phi.fun, old.phi = phi, 
            psi1 = exp(candidate), Nphys = Nphys, Ntune = Ntune)
        f.out <- p.eqn4.supp(D1 = D1, y = y, H1 = H1, lognormally.distributed = lognormally.distributed, 
            include.prior = include.prior, return.log = TRUE, 
            phi = phi.temp)
        if (do.print) {
            print(candidate)
            print(f.out)
        }
        if (do.filewrite) {
            save(phi.temp, f.out, ascii = TRUE, file = file.path(directory, 
                paste("stage1.", gsub("[ :]", "_", date()), sep = "")))
            save(phi.temp, f.out, ascii = TRUE, file = file.path(directory, 
                "stage1.latest"))
        }
        return(f.out)
    }
    psi1.start <- log(phi$psi1)
    e <- optim(psi1.start, f, method = method, control = list(fnscale = -1, 
        maxit = maxit, trace = 100))
    phi.out <- phi.change(phi.fun = phi.fun, old.phi = phi, psi1 = exp(e$par), Nphys = Nphys, Ntune = Ntune)
    return(invisible(phi.out))
}


## For the stage2 optimization we must define the H2 function, 
## which is the column vector of the regression function values 
## for the physical experiment data.  It is defined analogously
## to H1 above for constant mean models.

H2.ex<-function(D2)
{
out<-rep(1,dim(D2)[1])
dim(out)<-c(dim(D2)[1],1)
return(out)
}

## A second function that must be defined is the E.theta function.
## This computes the expectation of H1.ex with respect to theta.
## In our constant mean computer experiment H1.ex is a column vector 
## of all 1's so the expectation is H1.ex.  
## Note that this differs from the E.theta.toy function in 
## BACCO, which involved linear and other regression terms.


E.theta.ex <- function (D2 = NULL, H1 = NULL, x1 = NULL, x2 = NULL, phi, give.mean = TRUE) 
{
    if (give.mean) {
        m_theta <- phi$theta.apriori$mean
        return(H1(D1.fun(D2, t.vec = m_theta)))
    }
    else {
        out <- 0
        dim(out) = c(1,1)
        return(out)
    }
}


## Edash.theta.ex is similar to E.theta.ex, except it computes the 
## expectation of H1.ex with respect to Edash (see BACCO paper).
## In our constant mean computer experiment H1.ex is a column vector 
## of all 1's so the expectation is H1.ex.  
## Note that this differs from the E.theta.toy function in 
## BACCO, which involved linear and other regression terms.



Edash.theta.ex <- function (x, t.vec, k, H1, fast.but.opaque = FALSE, a = NULL, b = NULL, phi = NULL) 
{
    if (fast.but.opaque) {
        edash.mean <- a + crossprod(b, t.vec[k, ])
    }
    else {
        V_theta <- phi$theta.apriori$sigma
        m_theta <- phi$theta.apriori$mean
        omega_t <- phi$omega_t
        edash.mean <- solve(solve(V_theta) + 2 * omega_t, solve(V_theta, 
            m_theta) + 2 * crossprod(omega_t, t.vec[k, ]))
    }
    jj <- as.vector(edash.mean)
    names(jj) <- rownames(edash.mean)
    edash.mean <- jj
    return(H1(D1.fun(x, edash.mean)))
}


## Finally we must define the extractor function which splits 
## xcode.training into the values of the physical variables 
## and the values of the tuning variables.  Notice in our 
## example the values of the physical variable are in column 1
## and the values of the tuning variables in column 2.  In the 
## example below x.star contains the values of the physical 
## variable and t.vec the values of the tuning variables.

extractor.ex <- function(D1, Nphys, Ntune)
{
  Np1 = Nphys + 1
  Npt = Nphys + Ntune
  return(list(x.star=D1[, 1:Nphys, drop=FALSE],
              t.vec=D1[, Np1:Npt, drop=FALSE]))
}

## One of the functions stage2 calls is the function tt.fun
## When theta is 1-dim example, there is a slight modification 
## that we need to make.  

tt.fun <- function (D1, extractor, x.i, x.j, test.for.symmetry = FALSE, 
    method = 1, phi, Nphys, Ntune) 
{
    V_theta <- phi$theta.apriori$sigma
 
if (Ntune == 1) {
    dim(V_theta) <- c(1,1)
                }

    m_theta <- phi$theta.apriori$mean
    x.star <- extractor(D1, Nphys, Ntune)$x.star
    t.vec <- extractor(D1, Nphys, Ntune)$t.vec
    omega_x <- phi$omega_x
    omega_t <- phi$omega_t
    sigma1squared <- phi$sigma1squared
    A.lower <- phi$A.lower
    if (method == 0) {
        x.jframe <- kronecker(rep(1, nrow(t.vec)), t(x.j))
        x.iframe <- kronecker(rep(1, nrow(t.vec)), t(x.i))
        attributes(x.jframe) <- attributes(x.star)
        attributes(x.iframe) <- attributes(x.star)
        exp.x.jk <- exp(-dists.2frames(x.jframe, x.star, omega_x))
        exp.x.il <- exp(-dists.2frames(x.star, x.iframe, omega_x))
        exp.t.kl <- exp(-dists.2frames(t.vec, t.vec, omega_t)/2)
    }
    else {
        exp.x.jk <- kronecker(exp(-dists.2frames(t(x.j), x.star, 
            omega_x)), t(rep(1, nrow(D1))))
        exp.x.il <- kronecker(exp(-t(dists.2frames(x.star, t(x.i), 
            omega_x))), rep(1, nrow(D1)))
        exp.t.kl <- exp(-dists.2frames(t.vec, t.vec, omega_t)/2)
        rownames(exp.x.jk) <- rownames(D1)
        colnames(exp.x.jk) <- rownames(D1)
        rownames(exp.x.il) <- rownames(D1)
        colnames(exp.x.il) <- rownames(D1)
    }
    mtheta.minus.tmean <- matrix(0, nrow(x.star), nrow(x.star))
    for (k in 1:nrow(x.star)) {
        if (test.for.symmetry == TRUE) {
            l.vector <- 1:nrow(x.star)
        }
        else {
            l.vector <- k:nrow(x.star)
        }
        for (l in l.vector) {
            jj <- m_theta - (t.vec[k, ] + t.vec[l, ])/2
            A <- phi$A
            mtheta.minus.tmean[k, l] <- exp(-as.vector(crossprod(jj, 
                A) %*% jj)/2)
        }
    }
    if (test.for.symmetry == FALSE) {
        mtheta.minus.tmean <- symmetrize(mtheta.minus.tmean)
    }
    fish <- 1/sqrt(det(diag(nrow = nrow(V_theta)) + 4 * V_theta %*% 
        omega_t))
    out <- sigma1squared^2 * fish * exp.x.jk * exp.x.il * exp.t.kl * 
        mtheta.minus.tmean
    colnames(out) <- rownames(out)
    return(out)
}

## Next, we must modify the function ht.fun when theta is 
## 1-dimensional, just as we did tt.fun

ht.fun <- function (x.i, x.j, D1, extractor, Edash.theta, H1, fast.but.opaque = TRUE,  x.star = NULL, t.vec = NULL, phi, Nphys, Ntune) 
{
    V_theta <- phi$theta.apriori$sigma

if (Ntune == 1) {
    dim(V_theta) <- c(1,1)
                }

    m_theta <- phi$theta.apriori$mean
    omega_x <- phi$omega_x
    omega_t <- phi$omega_t
    if (fast.but.opaque == TRUE) {
        if (is.null(x.star) | is.null(t.vec)) {
            stop("argument 'fast.but.opaque' being TRUE means caller has to supply x.star and t.vec (use extractor.ex() for this)")
        }
    }
    else {
        x.star <- extractor(D1, Nphys, Ntune)$x.star
        t.vec <- extractor(D1, Nphys, Ntune)$t.vec
    }
    f1 <- function(k) {
        as.vector(quad.form(omega_x, x.i - x.star[k, ]))
    }
    matrix.in.f2 <- 2 * V_theta + solve(omega_t)
    f2 <- function(k) {
        as.vector(quad.form.inv(matrix.in.f2, m_theta - t.vec[k, 
            ]))
    }
    jj.a <- phi$a
    jj.b <- phi$b
    jj.c <- phi$c
    f3 <- function(k) {
        if (fast.but.opaque == TRUE) {
            return(Edash.theta(x = x.j, t.vec = t.vec, k = k, 
                H1 = H1, 
                phi = phi))
        }
        else {
            return(Edash.theta(x = x.j, t.vec = t.vec, k = k, 
                H1 = H1, phi = phi))
        }
    }
    out <- sapply(1:nrow(t.vec), function(i) {
        jj.c * exp(-f1(i) - f2(i)) * f3(i)
    })
    dim(out) = c(1,length(out))
    rownames(out) <- colnames(f3(1))
    colnames(out) <- rownames(t.vec)
return(t(out))
}



## Next we have to modify the function V.fun

V.fun <- function (D1, D2, H1, H2, extractor, E.theta, Edash.theta, give.answers = FALSE, test.for.symmetry = FALSE, phi, Nphys, Ntune) 
{
    lambda <- phi$lambda
    rho <- phi$rho
    x.star <- extractor(D1, Nphys, Ntune)$x.star
    t.vec <- extractor(D1, Nphys, Ntune)$t.vec
    v1.d1 <- V1(D1 = D1, other = NULL, phi = phi)
    v1.d1.inv <- solve(v1.d1)
    w1.d1 <- W1(D1 = D1, H1 = H1, phi = phi)
    h1.d1 <- H1(D1)
    w1.blah <- crossprod(w1.d1, crossprod(h1.d1, v1.d1.inv))
    v1.blah <- crossprod(v1.d1.inv, h1.d1) %*% w1.d1
    v1.blah.di.blah <- v1.blah %*% crossprod(h1.d1, v1.d1.inv)
    line1 <- line2 <- line3 <- line4 <- line5 <- line6 <- matrix(0, 
        nrow(D2), nrow(D2))
    tvec.arbitrary <- phi$theta.apriori$mean
    for (i in 1:nrow(D2)) {
        if (test.for.symmetry == FALSE) {
            j.vector <- i:nrow(D2)
        }
        else {
            j.vector <- 1:nrow(D2)
        }
        for (j in j.vector) {
            line1[i, j] <- V1(D1 = D1.fun(x.star = D2[i, ], t.vec = tvec.arbitrary),  other = D1.fun(x.star = D2[j, ], t.vec = tvec.arbitrary), 
phi = phi)
            jj.tt <- tt.fun(D1 = D1, extractor = extractor, x.i = D2[i, 
                ], x.j = D2[j, ], phi = phi, Nphys = Nphys, Ntune = Ntune)
            line2[i, j] <- tr(crossprod(v1.d1.inv, jj.tt))
            jj.hh <- hh.fun(x.i = D2[i, ], x.j = D2[j, ], E.theta = E.theta, 
                H1 = H1, phi = phi)
            line3[i, j] <- tr(crossprod(w1.d1, jj.hh))
            jj.th <- ht.fun(x.i = D2[j, ], x.j = D2[i, ], D1 = D1, 
                extractor = extractor, Edash.theta = Edash.theta, 
                H1 = H1, fast.but.opaque = TRUE, x.star = x.star, 
                t.vec = t.vec, phi = phi, Nphys = Nphys, Ntune = Ntune)
            line4[i, j] <- tr(w1.blah %*% jj.th)
            jj.ht <- ht.fun(x.i = D2[i, ], x.j = D2[j, ], D1 = D1, 
                extractor = extractor, Edash.theta = Edash.theta, 
                H1 = H1, fast.but.opaque = TRUE, x.star = x.star, 
                t.vec = t.vec, phi = phi, Nphys = Nphys, Ntune = Ntune)
            line5[i, j] <- tr(crossprod(v1.blah, jj.ht))
            line6[i, j] <- tr(crossprod(jj.tt, v1.blah.di.blah))
        }
    }
    C.m <- +line1 - line2 + line3 - line4 - line5 + line6
    if (test.for.symmetry == FALSE) {
        C.m <- symmetrize(C.m)
    }
    C.m <- rho^2 * C.m
    rownames(C.m) <- rownames(D2)
    colnames(C.m) <- rownames(D2)
    C.lambda <- lambda * diag(nrow = nrow(D2))
    C.v2 <- V2(D2, other = NULL, phi = phi)
    out <- C.lambda + C.v2 + C.m

    if (give.answers) {
        attributes(line1) <- attributes(out)
        attributes(line2) <- attributes(out)
        attributes(line3) <- attributes(out)
        attributes(line4) <- attributes(out)
        attributes(line5) <- attributes(out)
        attributes(line6) <- attributes(out)
        return(list(line1 = line1, line2 = line2, line3 = line3, 
            line4 = line4, line5 = line5, line6 = line6, v1.d1 = v1.d1, 
            v1.d1.inv = v1.d1.inv, w1.d1 = w1.d1, h1.d1 = h1.d1, 
            w1.blah = w1.blah, v1.blah = v1.blah, v1.blah.di.blah = v1.blah.di.blah, 
            C.m = C.m, C.lambda = C.lambda, C.v2 = C.v2, out = out))
    }
    else {
        return(out)
    }
}

## For one dimension we must modify the W2 function

W2 <- function (D2, H2, V, det = FALSE, Ntune) 
{
    out <- quad.form(solve(V), H2(D2))
if (length(out) == 1) {
                     dim(out) <-c(1,1)
                     }
if (Ntune == 1) {
                 dim(out) <-c(1,1)
                }
    if (det) {
        return(1/det(out))
    }
    else {
        return(solve(out))
    }
}



## Next, we need to modify the beta2hat.fun function for the 
## one-dimensional case

beta2hat.fun <- function (D1, D2, H1, H2, V = NULL, z, etahat.d2, extractor, E.theta, Edash.theta, phi, Nphys, Ntune) 
{
    rho <- phi$rho
    if (is.null(V)) {
        V <- V.fun(D1 = D1, D2 = D2, H1 = H1, H2 = H2, extractor = extractor, E.theta = E.theta, Edash.theta = Edash.theta, phi = phi, Nphys = Nphys, Ntune = Ntune)
    }
if (Ntune == 1) {
                 out <- crossprod(W2(D2, H2, V, Ntune = Ntune), 
                 crossprod(H2(D2),
                 solve(V, t(z - rho * etahat.d2))))
                 return(out)
                }
out <- crossprod(W2(D2, H2, V, Ntune = Ntune), crossprod(H2(D2), 
          solve(V, z - rho * etahat.d2)))
return(out)
}

## Next we modify the function t.fun

t.fun <- function (x, D1, extractor, phi, Nphys, Ntune) 
{
    sigma1squared <- phi$sigma1squared
    rho <- phi$rho
    omega_x <- phi$omega_x
    omega_t <- phi$omega_t
    x.star <- extractor(D1, Nphys, Ntune)$x.star
    t.vec <- extractor(D1, Nphys, Ntune)$t.vec
    m_theta <- phi$theta.apriori$mean
    V_theta <- phi$theta.apriori$sigma
    jj.mat <- solve(2 * V_theta + solve(omega_t))
    prod.1 <- sigma1squared/sqrt(det(diag(nrow = nrow(omega_t)) + 
        2 * crossprod(V_theta, omega_t)))
    out <- rep(NA, nrow(D1))
    for (j in 1:length(out)) {
        jj.x <- as.vector(x.star[j, ] - x)
        jj.m <- as.vector(m_theta - t.vec[j, ])
        prod.2 <- exp(-quad.form(omega_x, jj.x))
        prod.3 <- exp(-quad.form(jj.mat, jj.m))
        out[j] <- prod.1 * prod.2 * prod.3
    }
    names(out) <- rownames(t.vec)
    return(out)
}



## Next we modify the function etahat

etahat <- function (D1, D2, H1, y, E.theta, extractor, phi, Nphys, Ntune) 
{
    jj.hfun <- E.theta(D2 = D2, H1 = H1, phi = phi, give.mean = TRUE)
    jj.beta <- beta1hat.fun(D1 = D1, H1 = H1, y = y, phi = phi)
    bit1 <- drop(jj.hfun %*% jj.beta)
    b <- beta1hat.fun(D1 = D1, H1 = H1, y = y, phi = phi)
    f <- function(x) {
        t.fun(x, D1 = D1, extractor = extractor, phi = phi, Nphys = Nphys, Ntune = Ntune)
    }
    jj.tfun <- apply(D2, 1, f)
    bit2 <- crossprod(jj.tfun, solve(V1(D1, phi = phi), y - H1(D1) %*% 
        b))
    return(drop(bit1 + bit2))
}


## Next we modify the function p.page4

p.page4 <- function (D1, D2, H1, H2, V, y, z, E.theta, Edash.theta, extractor, include.prior = FALSE, lognormally.distributed = FALSE, return.log = FALSE, phi, Nphys, Ntune) 
{
    if (is.null(V)) {
        V <- V.fun(D1 = D1, D2 = D2, H1 = H1, H2 = H2, extractor = extractor, E.theta = E.theta, Edash.theta = Edash.theta, give.answers = FALSE, phi = phi, Nphys = Nphys, Ntune = Ntune)
    }
    etahat.d2 <- etahat(D1 = D1, D2 = D2, H1 = H1, y = y, E.theta = E.theta, extractor = extractor, phi = phi, Nphys, Ntune)
    beta2hat <- beta2hat.fun(D1 = D1, D2 = D2, H1 = H1, H2 = H2, 
        V = V, z = z, etahat.d2 = etahat.d2, extractor = extractor, 
        E.theta = E.theta, Edash.theta = Edash.theta, phi = phi, Nphys = Nphys, Ntune = Ntune)
    h2.d2 <- H2(D2)
    rho <- phi$rho
if (Ntune == 1) {
                 vec <-t(z) - h2.d2 %*% beta2hat - rho * etahat.d2
                }
else            {
                 vec <- z - h2.d2 %*% beta2hat - rho * etahat.d2
                } 

    if (return.log) {
        bit1 <- log(W2(D2, H2, V, det = TRUE, Ntune = Ntune))/2 - 
                       log(det(V))/2
        bit2 <- -0.5 * quad.form.inv(V, vec)
        if (include.prior) {
            return(bit1 + bit2 + log(prob.psi2(phi, lognormally.distributed = lognormally.distributed)))
        }
        else {
            return(bit1 + bit2)
        }
    }
    else {
        bit1 <- sqrt(W2(D2, H2, V, det = TRUE, Ntune = Ntune)/det(V))
        bit2 <- exp(-0.5 * quad.form.inv(V, vec))
        if (include.prior) {
            return(bit1 * bit2 * prob.psi2(phi, lognormally.distributed = lognormally.distributed))
        }
        else {
            return(bit1 * bit2)
        }
    }
}



## Now we must define the function stage2.  This is essentially 
## as in BACCO, but certain lines need to be changed to reflect 
## the dimensionality of our example.

stage2 <- function (D1, D2, H1, H2, y, z, maxit, method = "Nelder-Mead", 
    directory = ".", do.filewrite = FALSE, do.print = TRUE, extractor, 
    phi.fun, E.theta, Edash.theta, lognormally.distributed = FALSE, 
    use.standin = FALSE, phi, Nphys, Ntune) 
{
    if (do.filewrite & !isTRUE(file.info(directory)$isdir)) {
        stop("do.filewrite = TRUE; directory name supplied does not exist")
    }

Np3 = Nphys + 3

    f <- function(candidate) {
        phi.temp <- phi.change(phi.fun = phi.fun, old.phi = phi, 
            rho = exp(candidate[1]), lambda = exp(candidate[2]), 
            psi2 = exp(candidate[3:Np3]), Nphys = Nphys, Ntune = 
            Ntune)
        if (use.standin) {
            V.temp <- 1 + diag(3, nrow = nrow(D2))
        }
        else {
            V.temp <- V.fun(D1 = D1, D2 = D2, H1 = H1, H2 = H2, 
                extractor = extractor, E.theta = E.theta, Edash.theta = Edash.theta,  phi = phi.temp, Nphys = Nphys, Ntune = Ntune)
        }
        f.out <- drop(p.page4(D1 = D1, D2 = D2, H1 = H1, H2 = H2, 
            V = V.temp, z = z, y = y, E.theta = E.theta, Edash.theta = Edash.theta,  extractor = extractor, include.prior = FALSE, lognormally.distributed = lognormally.distributed,  return.log = TRUE, 
phi = phi.temp, Nphys = Nphys, Ntune = Ntune))
        if (do.print) {
            print(candidate)
            print(f.out)
        }
        if (do.filewrite) {
            save(phi.temp, f.out, ascii = TRUE, file = file.path(directory, 
                paste("stage2.", gsub("[ :]", "_", date()), sep = "")))
            save(phi.temp, f.out, ascii = TRUE, file = file.path(directory, 
                "stage2.latest"))
        }
        return(f.out)
    }
    rho.lambda.psi2.start <- log(c(phi$rho, phi$lambda, phi$psi2))
    e <- optim(rho.lambda.psi2.start, f, method = method, control = list(fnscale = -1, 
        trace = 100, maxit = maxit))
    jj <- exp(e$par)
    phi.out <- phi.change(phi.fun = phi.fun, old.phi = phi, rho = jj[1], 
        lambda = jj[2], psi2 = jj[3:Np3], Nphys = Nphys, Ntune = Ntune)
    return(invisible(phi.out))
}



## We next must modify the function stage 3.  Actually, we 
## need to modify the functions prob.theta, p.eqn8.supp, and 
## p.eqn8.supp.vector, which are called by stage 3.  The problem 
## is if theta is one-dimensional the function prob.theta 
## will return an error message.  prob.theta calls 
## the function dmvnorm, which uses the ncol function.
## Unfortunately, ncol requires a matrix or dataframe 
## as an object.  Likewise, p.eqn8.supp and 
## p.eqn8.supp.vector have problems when theta 
## is one-dimensional also.  When I tried to modify 
## these functions individually, I ran into difficulties 
## when I called stage3.  R seems to override my 
## redefinitions and use the original versions.
## To get around this, I have redefined prob.theta, 
## p.eqn8.supp, and p.eqn8.supp.vector inside 
## stage3.  This seems to avoid problems.  Note 
## that inside stage3 I renamed prob.theta p.theta.

stage3 <- function (D1, D2, H1, H2, d, maxit, trace = 100, method = "Nelder-Mead", 
    directory = ".", do.filewrite = FALSE, do.print = TRUE, include.prior = TRUE, 
    lognormally.distributed = FALSE, theta.start = NULL, phi) 
{

p.eqn8.supp <- function (theta, D1, D2, H1, H2, d, include.prior = FALSE, lognormally.distributed = FALSE, 
    return.log = FALSE, phi) 
{

p.eqn8.supp.vector <- function (theta, D1, D2, H1, H2, d, include.prior = FALSE, lognormally.distributed = FALSE, 
    return.log = FALSE, phi) 
{
p.theta <- function (theta, phi, lognormally.distributed = FALSE) 
{

if (length(theta) == 1) {
             return(dnorm(x = theta, mean = phi$theta.aprior$mean,
                    sd = sqrt(phi$theta.aprior$sigma)))
                        }
else   {
    if (lognormally.distributed) {
        if (any(x) < 0) {
            return(0)
                        }
        return(dmvnorm(x = log(theta), mean = log(phi$theta.aprior$mean), 
            sigma = phi$theta.aprior$sigma))
                                 }
       else {
        return(dmvnorm(x = theta, mean = phi$theta.aprior$mean, 
            sigma = phi$theta.aprior$sigma))
             }
          }
}


    betahat <- betahat.fun.koh(theta = theta, D1 = D1, D2 = D2, 
        H1 = H1, H2 = H2, phi = phi, d = d)
    Vd.theta <- Vd(theta = theta, D1 = D1, D2 = D2, phi = phi)
    det.W.theta <- W(theta = theta, D1 = D1, D2 = D2, H1 = H1, 
        H2 = H2, det = TRUE, phi = phi)
    H.theta <- H.fun(theta = theta, D1 = D1, D2 = D2, H1 = H1, 
        H2 = H2, phi = phi)
    mismatch <- d - H.theta %*% betahat
    if (return.log) {
        out <- -log(det(Vd.theta))/2 + log(det.W.theta)/2 - quad.form.inv(Vd.theta, 
            mismatch)/2
        if (include.prior) {
            out <- out + log(p.theta(theta, phi, lognormally.distributed = lognormally.distributed))
        }
    }
    else {
        out <- det(Vd.theta)^(-1/2) * det.W.theta^(1/2) * exp(-0.5 * 
            quad.form.inv(Vd.theta, mismatch))
        if (include.prior) {
            out <- out * p.theta(theta, phi, lognormally.distributed = lognormally.distributed)
        }
    }
    return(as.vector(out))
}


    if (is.vector(theta)) {
        return(p.eqn8.supp.vector(theta = theta, D1 = D1, D2 = D2, 
            H1 = H1, H2 = H2, d = d, include.prior = include.prior, 
            lognormally.distributed = lognormally.distributed, 
            return.log = return.log, phi = phi))
    }

else{
   return(0)
    }
}


    if (do.filewrite & !isTRUE(file.info(directory)$isdir)) {
        stop("do.filewrite = TRUE; directory name supplied does not exist")
    }
    f <- function(theta) {
        print(theta)
        f.out <- p.eqn8.supp(theta, D1 = D1, D2 = D2, H1 = H1, 
            H2 = H2, d = d, include.prior = include.prior, lognormally.distributed = lognormally.distributed, 
            return.log = TRUE, phi = phi)
        if (do.print) {
            print(theta)
            print(f.out)
        }
        if (do.filewrite) {
            save(theta, f.out, ascii = TRUE, file = file.path(directory, 
                paste("stage3.", gsub("[ :]", "_", date()), sep = "")))
            save(theta, f.out, ascii = TRUE, file = file.path(directory, 
                "stage3.latest"))
        }
        return(f.out)
    }
    if (is.null(theta.start)) {
        theta.start <- phi$theta.apriori$mean
    }
    e <- optim(theta.start, f, method = method, control = list(fnscale = -1, 
        maxit = maxit, trace = trace))
    optimal.theta <- e$par
    names(optimal.theta) <- names(phi$theta.apriori$mean)
    return(optimal.theta)
}