f <- function(x)
  {
    matrix(c(1, x[1], x[2], x[1]*x[1], x[1]*x[2], x[2]*x[2], x[1]^3), nrow=1)
  }

big.x <- function(doe)
  {
    t(apply(doe, 1, f))
  }

xtxinv <- function(big.x)
  {
    solve(t(big.x) %*% big.x)
  }

vary <- function(x, xtxinv)
  {
    (f(x) %*% xtxinv %*% t(f(x)))[1,1]
  }

iv <- function(m, doe)
  {
    sum(diag(m %*% xtxinv(big.x(doe))))
  }

x.pts <- seq(-1.0, 1.0, len=106)
plt.pts <- data.frame(x1=x.pts, x2=rep(x.pts, each=106))

n.terms <- 7
n.vars <- 2

m1 <- matrix(scan("doe2_moments.txt", sep="\t", what=double(0), quiet=TRUE),
             nrow=n.terms)

library(lattice)
# trellis.device(pdf, width=8, height=8, onefile=TRUE, file="doe2.pdf", version="1.4")
trellis.device(png, width=800, height=700, file="doe2-%02d.png")
for (n.runs in n.terms:25)
  {
    doe <- t(matrix(scan(paste("doe2/v", n.runs, "best", sep="."),
                         what=double(0), quiet=TRUE), 
                    nrow=n.vars))
    tmp.xtxinv <- xtxinv(big.x(doe))
    plt.pts$vary <- apply(plt.pts, 1, function(x) vary(x, tmp.xtxinv))
    main.title <- paste(n.runs, "Runs, IV =", 
                        formatC(iv(m1, doe), digits=9, format="f"))
    ## Now, make the plot
    plt <- levelplot(vary ~ x1*x2,
                     data=plt.pts,
                     contour=TRUE,
                     pretty=TRUE,
                     aspect=1,
                     xlab="x1",
                     ylab="x2",
                     main=main.title,
                     xlim=c(-1.05, 1.05),
                     ylim=c(-1.05, 1.05),
                     panel=function(x, y, subscripts, ...)
                     {
                       panel.levelplot(x, y, subscripts, ...)
                       doe <- round(doe*50) / 50
                       doe.u <- unique(doe)
                       doe.freq <-
                         apply(doe.u, 1,
                               function(x)
                               sum(apply(doe, 1, function(y) identical(x, y))))
                       panel.points(doe.u[,1], doe.u[,2], cex=3, pch=16, col="black")
                       panel.text(doe.u[,1], doe.u[,2], cex=1.2, labels=doe.freq,
                                  col="white")
                     })
    print(plt)
  }
dev.off()