t11d6 = read.table("table-11-6.txt", header = TRUE);
l = lm(y ~ x1 + x2 + I(x1^2) + I(x2^2) + x1:x2, t11d6);
print(summary(l));

ng = 31;
x1g = seq(from = min(t11d6$x1), to = max(t11d6$x1), length = ng);
x2g = seq(from = min(t11d6$x2), to = max(t11d6$x2), length = ng);
xi1g = seq(from = min(t11d6$xi1), to = max(t11d6$xi1), length = ng);
xi2g = seq(from = min(t11d6$xi2), to = max(t11d6$xi2), length = ng);

fit = matrix(0, ng, ng);
x1m = x1g[row(fit)];
x2m = x2g[col(fit)];

p = predict(l, data.frame(x1 = x1m[], x2 = x2m[]), se.fit = TRUE);
pdf("../response-fit-contour.pdf")
contour(xi1g, xi2g, matrix(p$fit, ng, ng),
        xlab = "Time", ylab = "Temperature", main = "Predicted Yield");
dev.off();

pdf("../response-fit-persp.pdf")
persp(xi1g, xi2g, matrix(p$fit, ng, ng), theta = -45, phi = 30,
      xlab = "Time", ylab = "Temperature", zlab = "Yield",
      expand = 0.6, r = 3, ticktype = "detailed");
dev.off();

pdf("../response-se-contour.pdf")
contour(xi1g, xi2g, matrix(p$se.fit, ng, ng),
        xlab = "Time", ylab = "Temperature", main = "Standard Error");
dev.off();

pdf("../response-se-persp.pdf")
persp(xi1g, xi2g, matrix(p$se.fit, ng, ng), theta = -45, phi = 30,
      xlab = "Time", ylab = "Temperature", zlab = "Standard Error",
      expand = 0.6, r = 3, ticktype = "detailed");
dev.off();

betaHat = coefficients(l);
b = betaHat[2:3];
B = matrix(c(betaHat[4], betaHat[6]/2, betaHat[6]/2, betaHat[5]), 2, 2);
cat("Stationary point:", -solve(B, b) / 2, "\n");

sl = summary(l)
V = sl$cov.unscaled;
msE = sl$sigma^2;
dfE = sl$df[2];

F = matrix(0, ng, ng);
for (i in 1:nrow(F)) {
        for (j in 1:ncol(F)) {
                L = rbind(
                        c(0, 1, 0, 2 * x1g[i],          0, x2g[j]),
                        c(0, 0, 1,          0, 2 * x2g[j], x1g[i])
                        );
                LbetaHat = L %*% betaHat;
                ssL = sum(LbetaHat * solve(L %*% V %*% t(L), LbetaHat));
                msL = ssL / 2;
                F[i, j] = msL / msE;
        }
}

pdf("../response-f-contour.pdf")
contour(xi1g, xi2g, F, col = "red",
        xlab = "Time", ylab = "Temperature", main = "F-Statistic");
contour(xi1g, xi2g, F, add = TRUE, lty = 2, levels = qf(.95, 2, dfE));
dev.off();

pdf("../response-f-persp.pdf")
persp(xi1g, xi2g, F, theta = -45, phi = 30,
      xlab = "Time", ylab = "Temperature", zlab = "F-Statistic",
      expand = 0.6, r = 3, ticktype = "detailed");
dev.off();

pdf("../response-f-detail-contour.pdf")
detail1 = xi1g > 85 & xi1g < 90;
detail2 = xi2g > 174 & xi2g < 180;
contour(xi1g[detail1], xi2g[detail2], F[detail1, detail2],
        col = "red",
        xlab = "Time", ylab = "Temperature", main = "F-Statistic Detail");
contour(xi1g, xi2g, F, add = TRUE, lty = 2, levels = qf(.95, 2, dfE));
dev.off();