---
title: "Calibration and Predictive Models"
author: "Ryan McClarren"
date: "November 18, 2014"
output: html_document
---

We're going to load in the shock breakout data from Stripling, McClarren, et al., ``A calibration and data assimilation method using the Bayesian MARS emulator''. Annals of Nuclear Energy, 52(C), 103–112.
```{r, include=FALSE}
require(ggplot2)
require(GGally)
require(tgp)
require(caret)
```
```{r}
fullData <- read.csv('CRASHBreakout.csv')
simulations <- fullData[,c("thickness","laser.energy","Be.gamma","wall.opacity","flux.limiter","breakout.time")]
measurements <- fullData[1:8,c("measure.1", "measure.2", "measure.3", "thickness", "laser.energy")]
```

The simulation has two physical parameters, $x$'s, Be disc thickness ($\mu$m) and laser energy (J), and 3 model parameters, $\theta$'s, Be $\gamma$ is a $\gamma$-law equation of state, wall opacity multiplier, and flux-limiter coefficient. There are `r length(simulations$thickness)` simulations in our data set. Looking at a pairs plot of the simulations, we can see that the breakout time in picoseconds (our quantity of interest) is strongly influenced by the disc thickness and the Be $\gamma$.
```{r, echo=FALSE}
ggpairs(simulations)
```

The first thing we will do is break the simulation data set into a training and test set. We'll use 80% of the data in the training set. This can be done simply using the <tt>createDataPartition</tt> function in the <tt>caret</tt> package:
```{r}
inTrain <- createDataPartition(y = simulations$breakout.time,times = 1, p = 0.8,list = F)
```
The next step is to use the <tt>bgp</tt> function (Bayesian Gaussian Process) from the <tt>tgp</tt> package.
```{r}
emulator <- bgp(X = simulations[inTrain,1:5], Z = simulations[inTrain,"breakout.time"], XX = simulations[-inTrain,1:5])
predDF <- data.frame(actual = simulations[inTrain,"breakout.time"], predicted = emulator$Zp.mean, q1 = emulator$Zp.q1, q2 = emulator$Zp.q2, type = "Training" )
predDF <- rbind(predDF, data.frame(actual = simulations[-inTrain,"breakout.time"], predicted = emulator$ZZ.mean, q1 = emulator$ZZ.q1, q2 = emulator$ZZ.q2, type = "Test" ))
```

Below is our comparison of the training and test data.
```{r}
ggplot(predDF,aes(x=predicted, y=actual, shape = type, color = type)) + geom_point() + geom_errorbarh(aes(xmin=predDF$q1, xmax = predDF$q2)) + geom_smooth(method="lm", se=F)
```

It seems like we can build a decent emulator for our simulations.  The root-mean square error for the training data is `r mean(sqrt((predDF$actual[ predDF$type =="Test"]-predDF$predicted[ predDF$type =="Test"])^2))` ps with out of an average breakout time of `r mean(predDF$actual[predDF$type =="Test"])`. Before we proceed, let's check density plots of the errors:
```{r}
ggplot(predDF,aes(x=actual-predicted, color=type)) + geom_density()
```

With this real data, there are some anomalies, but nothing invalidating the model. 

