text:
# Comparison of Nonparametric Bayesian Gaussian Process estimates to standard the Parametric Bayesian approach

```r
setwd("~/Documents/code/nonparametric-bayes/inst/examples/BUGS/")
```

Plotting and knitr options, (can generally be ignored)

```r
opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, cache=FALSE, comment=NA,
               fig.width=6, fig.height=4)

library(knitcitations)

library(ggplot2) # plotting
opts_knit$set(upload.fun = socialR::flickr.url)
theme_set(theme_bw(base_size=10))
theme_update(panel.background = element_rect(fill = "transparent", colour = NA),
             plot.background = element_rect(fill = "transparent", colour = NA))
cbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
               "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
```

Load necessary libraries,

```r
library(nonparametricbayes) # loads the rest as dependencies
```

### Model and parameters

Uses the model derived in `citet("10.1080/10236190412331335373")`, of a Ricker-like growth curve with an allee effect, defined in the pdgControl package,

```r
f <- RickerAllee
p <- c(2, 8, 5)
K <- p[2]
allee <- p[3]
```

Various parameters defining noise dynamics, grid, and policy costs.

```r
sigma_g <- 0.05
sigma_m <- 0.0
z_g <- function() rlnorm(1, 0, sigma_g)
z_m <- function() 1+(2*runif(1, 0,  1)-1) * sigma_m
x_grid <- seq(0, 1.5 * K, length=50)
h_grid <- x_grid
profit <- function(x,h) pmin(x, h)
delta <- 0.01
OptTime <- 70  # stationarity with unstable models is tricky thing
reward <- 0
xT <- 0
Xo <- K # observations start from
x0 <- Xo # simulation under policy starts from
Tobs <- 20
```

### Sample Data

```r
  set.seed(1234)
  #harvest <- sort(rep(seq(0, .5, length=7), 5))
  x <- numeric(Tobs)
  x[1] <- 6
  nz <- 1
  for(t in 1:(Tobs-1))
    x[t+1] = z_g() * f(x[t], h=0, p=p)
  obs <- data.frame(x = c(rep(0,nz),
                          pmax(rep(0,Tobs-1), x[1:(Tobs-1)])),
                    y = c(rep(0,nz),
                          x[2:Tobs]))
```

## Maximum Likelihood

```r
alt <- par_est(obs,  init = c(r = p[1],
                              K = mean(obs$x[obs$x>0]),
                              s = sigma_g))
est <- par_est_allee(obs, f, p,
                     init = c(r = p[1] + 1,
                              K = p[2] + 2,
                              C = p[3] + 2,
                              s = sigma_g))
```

## Non-parametric Bayes

```r
#inv gamma has mean b / (a - 1) (assuming a>1) and variance b ^ 2 / ((a - 2) * (a - 1) ^ 2) (assuming a>2)
s2.p <- c(5,5)
d.p = c(10, 1/0.1)
```

Estimate the Gaussian Process (nonparametric Bayesian fit)

```r
gp <- gp_mcmc(obs$x, y=obs$y, n=1e5, s2.p = s2.p, d.p = d.p)
```

![plot of chunk unnamed-chunk-2](http://farm8.staticflickr.com/7322/8717734359_7715ddc40e_o.png)

```r
gp_dat <- gp_predict(gp, x_grid, burnin=1e4, thin=300)
```

Show traces and posteriors against priors

```r
plots <- summary_gp_mcmc(gp)
```

![plot of chunk gp_traces_densities](http://farm8.staticflickr.com/7437/8718854190_d70491a232_o.png) ![plot of chunk gp_traces_densities](http://farm8.staticflickr.com/7435/8717735439_7fab806a68_o.png)

```r
# Summarize the GP model
tgp_dat <-
    data.frame(  x = x_grid,
                 y = gp_dat$E_Ef,
                 ymin = gp_dat$E_Ef - 2 * sqrt(gp_dat$E_Vf),
                 ymax = gp_dat$E_Ef + 2 * sqrt(gp_dat$E_Vf) )
```

### Parametric Bayes

We initiate the MCMC chain (`init_p`) using the true values of the parameters `p` from the simulation.  While impossible in real data, this gives the parametric Bayesian approach the best chance at succeeding.  `y` is the timeseries (recall `obs` has the $x_t$, $x_{t+1}$ pairs)

```r
# a bit unfair to start with the correct values, but anyhow...
init_p = p
names(init_p) = c("r0", "K", "theta")
y <- obs$x[-1]
N=length(y);
```

We'll be using the JAGS Gibbs sampler, a recent open source BUGS implementation with an R interface that works on most platforms.  We initialize the usual MCMC parameters; see `?jags` for details.

```r
jags.data <- list("N","y")
n.chains = 1
n.iter = 40000
n.burnin = floor(10000)
n.thin = max(1, floor(n.chains * (n.iter - n.burnin)/1000))
```

The actual model is defined in a `model.file` that contains an R function that is automatically translated into BUGS code by *R2WinBUGS*.  The file defines the priors and the model, as seen when read in here

```r
cat(readLines(con="bugmodel-UPrior.txt"), sep="\n")
```

```
model{

K     ~ dunif(0.01, 40.0)
logr0    ~ dunif(-6.0, 6.0)
logtheta ~ dunif(-6.0, 6.0)
stdQ ~ dunif(0.0001,100)
stdR ~ dunif(0.0001,100)
# JAGS notation, mean, and precision ( reciprical of the variance, 1/sigma^2)
iQ <- 1/(stdQ*stdQ);
iR <- 1/(stdR*stdR);

r0 <- exp(logr0)
theta <- exp(logtheta)

x[1] ~ dunif(0,10)

for(t in 1:(N-1)){
  mu[t] <- x[t] * exp(r0 * (1 - x[t]/K)* (x[t] - theta) / K )
  x[t+1] ~ dnorm(mu[t],iQ)
}

for(t in 1:(N)){
  y[t] ~ dnorm(x[t],iR)
}

}
```

We define which parameters to keep track of, and set the initial values of parameters in the transformed space used by the MCMC.  We use logarithms to maintain strictly positive values of parameters where appropriate.

```r
# Uniform priors on standard deviation terms
jags.params=c("K","logr0","logtheta","stdQ", "stdR")
jags.inits <- function(){
  list("K"=init_p["K"],"logr0"=log(init_p["r0"]),"logtheta"=log(init_p["theta"]), "stdQ"=sqrt(0.05),"stdR"=sqrt(0.1),"x"=y,.RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}

set.seed(12345)

time_jags <- system.time(
  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params, n.chains=n.chains,
                  n.iter=n.iter, n.thin=n.thin, n.burnin=n.burnin,model.file="bugmodel-UPrior.txt")
)
```

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 185

Initializing model

  |
  |                                                  |   0%
  |
  |++++                                              |   8%
  |
  |++++++++                                          |  16%
  |
  |++++++++++++                                      |  24%
  |
  |++++++++++++++++                                  |  32%
  |
  |++++++++++++++++++++                              |  40%
  |
  |++++++++++++++++++++++++                          |  48%
  |
  |++++++++++++++++++++++++++++                      |  56%
  |
  |++++++++++++++++++++++++++++++++                  |  64%
  |
  |++++++++++++++++++++++++++++++++++++              |  72%
  |
  |++++++++++++++++++++++++++++++++++++++++          |  80%
  |
  |++++++++++++++++++++++++++++++++++++++++++++      |  88%
  |
  |++++++++++++++++++++++++++++++++++++++++++++++++  |  96%
  |
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%

  |
  |                                                  |   0%
  |
  |*                                                 |   3%
  |
  |***                                               |   5%
  |
  |****                                              |   8%
  |
  |*****                                             |  11%
  |
  |*******                                           |  13%
  |
  |********                                          |  16%
  |
  |*********                                         |  19%
  |
  |***********                                       |  21%
  |
  |************                                      |  24%
  |
  |*************                                     |  27%
  |
  |***************                                   |  29%
  |
  |****************                                  |  32%
  |
  |*****************                                 |  35%
  |
  |*******************                               |  37%
  |
  |********************                              |  40%
  |
  |*********************                             |  43%
  |
  |***********************                           |  45%
  |
  |************************                          |  48%
  |
  |*************************                         |  51%
  |
  |***************************                       |  53%
  |
  |****************************                      |  56%
  |
  |*****************************                     |  59%
  |
  |*******************************                   |  61%
  |
  |********************************                  |  64%
  |
  |*********************************                 |  67%
  |
  |***********************************               |  69%
  |
  |************************************              |  72%
  |
  |*************************************             |  75%
  |
  |***************************************           |  77%
  |
  |****************************************          |  80%
  |
  |*****************************************         |  83%
  |
  |*******************************************       |  85%
  |
  |********************************************      |  88%
  |
  |*********************************************     |  91%
  |
  |***********************************************   |  93%
  |
  |************************************************  |  96%
  |
  |************************************************* |  99%
  |
  |**************************************************| 100%
```

```r
time_jags <- unname(time_jags["elapsed"]);
```

#### Convergence diagnostics for parametric bayes

```r
jags_matrix <- as.data.frame(as.mcmc.bugs(jagsfit$BUGSoutput))
par_posteriors <- melt(cbind(index = 1:dim(jags_matrix)[1], jags_matrix), id = "index")
```

```r
# Traces
ggplot(par_posteriors) + geom_line(aes(index, value)) + facet_wrap(~ variable, scale="free", ncol=1)
```

![plot of chunk parametric_bayes_traces](http://farm8.staticflickr.com/7446/8718855418_3973517549_o.png)

```r
## priors (untransformed variables)
K_prior <- function(x) dunif(x, 0.01, 40)
logr_prior <- function(x) dunif(x, -6, 6)
logtheta_prior <- function(x) dunif(x, -6, 6)
stdQ_prior <- function(x) dunif(x, 0.001, 100)
stdR_prior <- function(x) dunif(x, 0.001, 100)

par_priors <- list(K = K_prior, deviance = function(x) 0 * x,
                   logr0 = logr_prior, logtheta = logtheta_prior,
                   stdQ = stdQ_prior, stdR = stdR_prior)

par_prior_curves <- ddply(par_posteriors, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})

# posterior distributions
ggplot(par_posteriors, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
  geom_line(data=par_prior_curves, aes(x=value, y=density), col="red") +
  facet_wrap(~ variable, scale="free", ncol=2)
```

![plot of chunk parametric_bayes_posteriors](http://farm8.staticflickr.com/7323/8718855656_60293c33b3_o.png)

```r
# um, cleaner if we were just be using the long form, par_posterior
mcmc <- as.mcmc(jagsfit)
mcmcall <- mcmc[,-2]
who <- colnames(mcmcall)
who
```

```
[1] "K"        "logr0"    "logtheta" "stdQ"     "stdR"
```

```r
mcmcall <- cbind(mcmcall[,1],mcmcall[,2],mcmcall[,3],mcmcall[,4],mcmcall[,5])
colnames(mcmcall) <- who

pardist <- mcmcall
pardist[,2] = exp(pardist[,2]) # transform model parameters back first
pardist[,3] = exp(pardist[,3])

bayes_coef <- apply(pardist,2,mean)
bayes_pars <- unname(c(bayes_coef[2], bayes_coef[1], bayes_coef[3]))
bayes_pars
```

```
[1]  0.1662 17.8364  2.9134
```

```r
par_bayes_means <- sapply(x_grid, f, 0, bayes_pars)
```

## Parametric Bayes based on the structurally wrong model

We initiate the MCMC chain (`init_p`) using the true values of the parameters `p` from the simulation.  While impossible in real data, this gives the parametric Bayesian approach the best chance at succeeding.  `y` is the timeseries (recall `obs` has the $x_t$, $x_{t+1}$ pairs)

```r
init_p = p
names(init_p) = c("r0", "K")
y <- obs$x[-1]
N=length(y);
```

We'll be using the JAGS Gibbs sampler, a recent open source BUGS implementation with an R interface that works on most platforms.  We initialize the usual MCMC parameters; see `?jags` for details.

```r
jags.data <- list("N","y")
n.chains = 1
n.iter = 400000
n.burnin = floor(10000)
n.thin = max(1, floor(n.chains * (n.iter - n.burnin)/1000))
```

The actual model is defined in a `model.file` that contains an R function that is automatically translated into BUGS code by *R2WinBUGS*.  The file defines the priors and the model, as seen when read in here

```r
cat(readLines(con="ricker-UPrior.txt"), sep="\n")
```

```
model{

K     ~ dunif(0.01, 40.0)
logr0    ~ dunif(-6.0, 6.0)
stdQ ~ dunif(0.0001,100)
stdR ~ dunif(0.0001,100)
# JAGS notation, mean, and precision ( reciprical of the variance, 1/sigma^2)
iQ <- 1/(stdQ*stdQ);
iR <- 1/(stdR*stdR);

r0 <- exp(logr0)

x[1] ~ dunif(0,10)

for(t in 1:(N-1)){
  mu[t] <- x[t] * exp(r0 * (1 - x[t]/K) / K )
  x[t+1] ~ dnorm(mu[t],iQ)
}

for(t in 1:(N)){
  y[t] ~ dnorm(x[t],iR)
}

}
```

We define which parameters to keep track of, and set the initial values of parameters in the transformed space used by the MCMC.  We use logarithms to maintain strictly positive values of parameters where appropriate.

```r
# Uniform priors on standard deviation terms
jags.params=c("K","logr0", "stdQ", "stdR")
jags.inits <- function(){
  list("K"=init_p["K"],"logr0"=log(init_p["r0"]), "stdQ"=sqrt(0.05),"stdR"=sqrt(0.1),"x"=y,.RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}

set.seed(12345)

time_jags <- system.time(
  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params, n.chains=n.chains,
                  n.iter=n.iter, n.thin=n.thin, n.burnin=n.burnin,model.file="ricker-UPrior.txt")
)
```

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 165

Initializing model

  |
  |                                                  |   0%
  |
  |++++++++++++++++++++++++++++++++++++++++          |  80%
  |
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%

  |
  |                                                  |   0%
  |
  |*                                                 |   2%
  |
  |**                                                |   4%
  |
  |***                                               |   6%
  |
  |****                                              |   8%
  |
  |*****                                             |  10%
  |
  |******                                            |  12%
  |
  |*******                                           |  14%
  |
  |********                                          |  16%
  |
  |*********                                         |  18%
  |
  |**********                                        |  21%
  |
  |***********                                       |  23%
  |
  |************                                      |  25%
  |
  |*************                                     |  27%
  |
  |**************                                    |  29%
  |
  |***************                                   |  31%
  |
  |****************                                  |  33%
  |
  |*****************                                 |  35%
  |
  |******************                                |  37%
  |
  |*******************                               |  39%
  |
  |*********************                             |  41%
  |
  |**********************                            |  43%
  |
  |***********************                           |  45%
  |
  |************************                          |  47%
  |
  |*************************                         |  49%
  |
  |**************************                        |  51%
  |
  |***************************                       |  53%
  |
  |****************************                      |  55%
  |
  |*****************************                     |  57%
  |
  |******************************                    |  59%
  |
  |*******************************                   |  62%
  |
  |********************************                  |  64%
  |
  |*********************************                 |  66%
  |
  |**********************************                |  68%
  |
  |***********************************               |  70%
  |
  |************************************              |  72%
  |
  |*************************************             |  74%
  |
  |**************************************            |  76%
  |
  |***************************************           |  78%
  |
  |****************************************          |  80%
  |
  |*****************************************         |  82%
  |
  |******************************************        |  84%
  |
  |*******************************************       |  86%
  |
  |********************************************      |  88%
  |
  |*********************************************     |  90%
  |
  |**********************************************    |  92%
  |
  |***********************************************   |  94%
  |
  |************************************************  |  96%
  |
  |************************************************* |  98%
  |
  |**************************************************| 100%
```

```r
time_jags <- unname(time_jags["elapsed"]);
```

#### Convergence diagnostics for parametric bayes Ricker model

```r
jags_matrix <- as.data.frame(as.mcmc.bugs(jagsfit$BUGSoutput))
par_posteriors <- melt(cbind(index = 1:dim(jags_matrix)[1], jags_matrix), id = "index")

# Traces
ggplot(par_posteriors) + geom_line(aes(index, value)) +
  facet_wrap(~ variable, scale="free", ncol=1)
```

![plot of chunk ricker_traces](http://farm8.staticflickr.com/7399/8718861954_b710a74a09_o.png)

```r
## priors (untransformed variables)
K_prior <- function(x) dunif(x, 0.01, 40)
logr_prior <- function(x) dunif(x, -6, 6)
stdQ_prior <- function(x) dunif(x, 0.001, 100)
stdR_prior <- function(x) dunif(x, 0.001, 100)

par_priors <- list(K = K_prior, deviance = function(x) 0 * x, logr0 = logr_prior, stdQ = stdQ_prior, stdR = stdR_prior)

par_prior_curves <- ddply(par_posteriors, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})

# posterior distributions
ggplot(par_posteriors, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
  geom_line(data=par_prior_curves, aes(x=value, y=density), col="red") +
  facet_wrap(~ variable, scale="free", ncol=2)
```

![plot of chunk ricker_posteriors](http://farm8.staticflickr.com/7403/8717743447_04d748d91c_o.png)

```r
ricker_pardist <- jags_matrix[! names(jags_matrix) == "deviance" ]
ricker_pardist[,"logr0"] = exp(ricker_pardist[,"logr0"]) # transform model parameters back first

posterior.mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

apply(ricker_pardist,2,mean)
```

```
      K   logr0    stdQ    stdR
19.2188  0.5702  0.3379  0.1492
```

```r
bayes_coef <- apply(ricker_pardist,2, posterior.mode) # much better estimates
ricker_bayes_pars <- unname(c(bayes_coef[2], bayes_coef[1]))
ricker_bayes_pars
```

```
[1] 0.7731 8.2288
```

## Write the external bugs file

```r
logr0_prior_p <- c(-6.0, 6.0)
logtheta_prior_p <- c(-6.0, 6.0)
logK_prior_p <- c(-6.0, 6.0)
stdQ_prior_p <- c(0.0001, 100)
stdR_prior_p <- c(0.0001, 100)

bugs.model <-
paste(sprintf(
"model{
  logr0    ~ dunif(%s, %s)
  logtheta    ~ dunif(%s, %s)
  logK    ~ dunif(%s, %s)
  stdQ ~ dunif(%s, %s)
  stdR ~ dunif(%s, %s)",
  logr0_prior_p[1], logr0_prior_p[2],
  logtheta_prior_p[1], logtheta_prior_p[2],
  logK_prior_p[1], logK_prior_p[2],
  stdQ_prior_p[1], stdQ_prior_p[2],
  stdR_prior_p[1], stdR_prior_p[2]),
  "

  iQ <- 1 / (stdQ * stdQ);
  iR <- 1 / (stdR * stdR);

  r0 <- exp(logr0)
  theta <- exp(logtheta)
  K <- exp(logK)

  x[1] ~ dunif(0, 10)

  for(t in 1:(N-1)){
    mu[t] <- r0 * pow(abs(x[t]), theta) / (1 + pow(abs(x[t]), theta) / K)
    x[t+1] ~ dnorm(mu[t], iQ)
  }

  for(t in 1:(N)){
    y[t] ~ dnorm(x[t], iR)
  }
}")

writeLines(bugs.model, "myers.bugs")
```

```r

## priors (untransformed variables)
logK_prior     <- function(x) dunif(x, logK_prior_p[1], logK_prior_p[2])
logr_prior     <- function(x) dunif(x, logr0_prior_p[1], logr0_prior_p[2])
logtheta_prior <- function(x) dunif(x, logtheta_prior_p[1], logtheta_prior_p[2])
stdQ_prior     <- function(x) dunif(x, stdQ_prior_p[1], stdQ_prior_p[2])
stdR_prior     <- function(x) dunif(x, stdR_prior_p[1], stdR_prior_p[2])
```

### Parametric Bayes based on myers model

Contains the

```r

init_p = c(1, 1, 1)
names(init_p) = c("r0", "theta", "K")
y <- obs$x[-1]
N=length(y);
```

We'll be using the JAGS Gibbs sampler, a recent open source BUGS implementation with an R interface that works on most platforms.  We initialize the usual MCMC parameters; see `?jags` for details.

```r
jags.data <- list("N","y")
n.chains = 1
n.iter = 40000
n.burnin = floor(10000)
n.thin = max(1, floor(n.chains * (n.iter - n.burnin)/1000))
```

```r
jags.params=c("logr0", "logtheta", "logK", "stdQ", "stdR")
jags.inits <- function(){
  list("logr0"=log(init_p["r0"]), "logtheta" = log(init_p["theta"]),
       "logK"=log(init_p["K"]), "stdQ"=sqrt(0.05),
       "stdR"=sqrt(0.1),"x"=y,.RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
```

```r
set.seed(12345)
time_jags <- system.time(
  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
                  n.chains=n.chains, n.iter=n.iter, n.thin=n.thin,
                  n.burnin=n.burnin, model.file="myers.bugs")
)
```

```
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 166

Initializing model

  |
  |                                                  |   0%
Deleting model
```

```
Error: Error in node x[15] Failure to calculate log density
```

```
Timing stopped at: 0.024 0 0.024
```

```r
time_jags <- unname(time_jags["elapsed"]);
```

Convergence diagnostics for parametric bayes

```r
jags_matrix <- as.data.frame(as.mcmc.bugs(jagsfit$BUGSoutput))
par_posteriors <- melt(cbind(index = 1:dim(jags_matrix)[1], jags_matrix), id = "index")

## Traces
ggplot(par_posteriors) + geom_line(aes(index, value)) +
  facet_wrap(~ variable, scale="free", ncol=1)
```

![plot of chunk unnamed-chunk-22](http://farm8.staticflickr.com/7298/8718862826_99d76f0a19_o.png)

```r
## prior curve functions, must be in order
par_priors <- list( deviance = function(x) 0 * x, logK = logK_prior,
                    logr0 = logr_prior, logtheta = logtheta_prior,
                    stdQ = stdQ_prior, stdR = stdR_prior)

## calculate points along the prior curve
par_prior_curves <- ddply(par_posteriors, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = par_priors[[dd$variable[1]]](grid))
})
```

```r
## posterior distributions
ggplot(par_posteriors, aes(value)) +
  stat_density(geom="path", position="identity", alpha=0.7) +
  geom_line(data=par_prior_curves, aes(x=value, y=density), col="red") +
  facet_wrap(~ variable, scale="free", ncol=2)
```

![plot of chunk unnamed-chunk-24](http://farm8.staticflickr.com/7305/8718864428_d2183d83ff_o.png)

```r
myers_pardist <- as.matrix(jags_matrix[2:6])
```

```
Error: undefined columns selected
```

```r
myers_pardist[,1] = exp(myers_pardist[,1]) # transform model parameters back first
myers_pardist[,2] = exp(myers_pardist[,2]) # transform model parameters back first
myers_pardist[,3] = exp(myers_pardist[,3]) # transform model parameters back first
colnames(myers_pardist) = c("K", "r0", "theta", "stdQ", "stdR")

apply(myers_pardist,2,mean)
```

```
         K         r0      theta       stdQ       stdR
1.310e+172  2.243e+01  1.834e+00  3.232e-01  1.778e-01
```

```r
bayes_coef <- apply(myers_pardist,2, posterior.mode) # much better estimates

myers_bayes_pars <- unname(c(bayes_coef[2], bayes_coef[3], bayes_coef[1]))
```

### Phase-space diagram of the expected dynamics

```r
true_means <- sapply(x_grid, f, 0, p)
alt_means <- sapply(x_grid, alt$f, 0, ricker_bayes_pars[c(1,2)])
est_means <- sapply(x_grid, est$f, 0, est$p)
par_bayes_means <- sapply(x_grid, f, 0, bayes_pars)

myers_means <- sapply(x_grid, Myer_harvest, 0, myers_bayes_pars)

models <- data.frame(x=x_grid, GP=tgp_dat$y, True=true_means,
                     MLE=est_means, Ricker=alt_means,
                     Parametric.Bayes = par_bayes_means,
                     Myers = myers_means)

models <- melt(models, id="x")
names(models) <- c("x", "method", "value")
```

```r
plot_gp <- ggplot(tgp_dat) + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") +
    geom_line(data=models, aes(x, value, col=method), lwd=1, alpha=0.8) +
    geom_point(data=obs, aes(x,y), alpha=0.8) +
    xlab(expression(X[t])) + ylab(expression(X[t+1])) +
    scale_colour_manual(values=cbPalette)
print(plot_gp)
```

![plot of chunk Figure1](http://farm8.staticflickr.com/7460/8718864688_a9b1f01827_o.png)

## Optimal policies by value iteration

Compute the optimal policy under each model using stochastic dynamic programming. We begin with the policy based on the GP model,

```r
MaxT = 1000
# uses expected values from GP, instead of integrating over posterior
#matrices_gp <- gp_transition_matrix(gp_dat$E_Ef, gp_dat$E_Vf, x_grid, h_grid)

# Integrate over posteriors
matrices_gp <- gp_transition_matrix(gp_dat$Ef_posterior, gp_dat$Vf_posterior, x_grid, h_grid)

# Solve the SDP using the GP-derived transition matrix
opt_gp <- value_iteration(matrices_gp, x_grid, h_grid, MaxT, xT, profit, delta, reward)
```

Determine the optimal policy based on the true and MLE models

```r
matrices_true <- f_transition_matrix(f, p, x_grid, h_grid, sigma_g)
opt_true <- value_iteration(matrices_true, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)

matrices_estimated <- f_transition_matrix(est$f, est$p, x_grid, h_grid, est$sigma_g)
opt_estimated <- value_iteration(matrices_estimated, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
```

Determine the optimal policy based on parametric Bayesian model

```r
allen_f <- function(x,h,p) unname(f(x,h,p[c(2, 1, 3)]))
matrices_par_bayes <- parameter_uncertainty_SDP(allen_f, x_grid, h_grid, pardist, 4)
opt_par_bayes <- value_iteration(matrices_par_bayes, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
```

Bayesian Ricker

```r
ricker_f <- function(x, h, p) alt$f(x, h, unname(p[c(2, 1)]))
matrices_alt <- parameter_uncertainty_SDP(ricker_f, x_grid, h_grid, as.matrix(ricker_pardist), 3)
```

![plot of chunk unnamed-chunk-30](http://farm8.staticflickr.com/7389/8718879588_3495523ab2_o.png)

```r
opt_alt <- value_iteration(matrices_alt, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
```

Bayesian Myers model

```r
myers_f <- function(x,h,p) Myer_harvest(x, h, p[c(2, 3, 1)])
matrices_myers <- parameter_uncertainty_SDP(myers_f, x_grid, h_grid, as.matrix(myers_pardist), 4)
myers_alt <- value_iteration(matrices_myers, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
```

Assemble the data

```r
OPT = data.frame(GP = opt_gp$D, True = opt_true$D, MLE = opt_estimated$D, Ricker = opt_alt$D, Parametric.Bayes = opt_par_bayes$D, Myers = myers_alt$D)
colorkey=cbPalette
names(colorkey) = names(OPT)
```

## Graph of the optimal policies

```r
policies <- melt(data.frame(stock=x_grid, sapply(OPT, function(x) x_grid[x])), id="stock")
names(policies) <- c("stock", "method", "value")

ggplot(policies, aes(stock, stock - value, color=method)) +
  geom_line(lwd=1.2, alpha=0.8) + xlab("stock size") + ylab("escapement")  +
  scale_colour_manual(values=colorkey)
```

![plot of chunk Figure2](http://farm8.staticflickr.com/7365/8717765091_59e0191a22_o.png)

## Simulate 100 realizations managed under each of the policies

```r
sims <- lapply(OPT, function(D){
  set.seed(1)
  lapply(1:100, function(i)
    ForwardSimulate(f, p, x_grid, h_grid, x0, D, z_g, profit=profit, OptTime=OptTime)
  )
})

dat <- melt(sims, id=names(sims[[1]][[1]]))
dt <- data.table(dat)
setnames(dt, c("L1", "L2"), c("method", "reps"))
# Legend in original ordering please, not alphabetical:
dt$method = factor(dt$method, ordered=TRUE, levels=names(OPT))
```

```r
ggplot(dt) +
  geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1) +
  scale_colour_manual(values=colorkey, guide = guide_legend(override.aes = list(alpha = 1)))
```

![plot of chunk Figure3](http://farm8.staticflickr.com/7398/8717766321_60e5ddc134_o.png)

```r
Profit <- dt[, sum(profit), by=c("reps", "method")]
Profit[, mean(V1), by="method"]
```

```
             method     V1
1:               GP 34.082
2:             True 34.256
3:              MLE  8.000
4:           Ricker  6.857
5: Parametric.Bayes 20.589
6:            Myers  8.000
```

```r
ggplot(Profit, aes(V1)) + geom_histogram() +
  facet_wrap(~method, scales = "free_y") + guides(legend.position = "none")
```

![plot of chunk totalprofits](http://farm8.staticflickr.com/7289/8717766899_5361fa4111_o.png)


text:
# Comparison of Nonparametric Bayesian Gaussian Process estimates to standard the Parametric Bayesian approach

Plotting and knitr options, (can generally be ignored)

```{r plotting-options, cache=FALSE, message=FALSE, warning=FALSE, include=FALSE}
source("~/.knitr_defaults.R")
#opts_knit$set(upload.fun = socialR::flickr.url)
library(knitcitations)
library(nonparametricbayes)
opts_chunk$set(external=TRUE)
read_chunk("external-chunks.R")
```

```{r posterior-mode}
```

### Model and parameters

```{r stateeq-local}
f <- Myers
p <- c(1.5, 2.5, 8)
K <- 5  # approx, a li'l' less
allee <- 1.2 # approx, a li'l' less
```

Various parameters defining noise dynamics, grid, and policy costs.

```{r sdp-pars, dependson="stateeq-local"}
```

### Sample Data

```{r obs, dependson="sdp-pars"}
```

## Maximum Likelihood

```{r mle, dependson="obs"}
```

Mean predictions

```{r mle-output, dependson="mle"}
```

## Non-parametric Bayes

```{r gp-priors}
```

Estimate the Gaussian Process (nonparametric Bayesian fit)

```{r gp, dependson=c("gp-priors", "obs")}
```

Show traces and posteriors against priors

```{r gp_traces_densities, dependson="gp"}
```

```{r gp-output, dependson="gp"}
```

Parametric Bayesian Models
---------------------------

We use the JAGS Gibbs sampler, a recent open source BUGS
implementation with an R interface that works on most platforms.
We initialize the usual MCMC parameters; see `?jags` for details.

All parametric Bayesian estimates use the following basic parameters for the JAGS MCMC:

```{r jags-setup}
```

We will use the same priors for process and observation noise in each model,

```{r common-priors}
```

### Parametric Bayes of correct (Allen) model

We initiate the MCMC chain (`init_p`) using the true values of the
parameters `p` from the simulation.  While impossible in real data, this
gives the parametric Bayesian approach the best chance at succeeding.
`y` is the timeseries (recall `obs` has the $x_t$, $x_{t+1}$ pairs)

The actual model is defined in a `model.file` that contains an R function
that is automatically translated into BUGS code by *R2WinBUGS*.  The file
defines the priors and the model. We write the file from R as follows:

```{r allen-model}
```

Write the priors into a list for later reference

```{r allen-priors, dependson="common-priors"}
```

We define which parameters to keep track of, and set the initial values of
parameters in the transformed space used by the MCMC.  We use logarithms
to maintain strictly positive values of parameters where appropriate.

```{r allen-mcmc, dependson=c("allen-model", "allen-pars", "jags-setup"), results="hide"}
```

#### Convergence diagnostics for Allen model

R notes: this strips classes from the `mcmc.list` object (so that we have list of matrices; objects that `reshape2::melt` can handle intelligently), and then combines chains into one array. In this array each parameter is given its value at each sample from the posterior (index) for each chain.

```{r allen-traces, dependson="allen-mcmc"}
```

```{r allen-posteriors, dependson=c("allen-traces", "allen-priors")}
```

Reshape the posterior parameter distribution data, transform back into original space, and calculate the mean parameters and mean function

```{r allen-output, dependson=c("posterior-mode", "allen-traces")}
```

## Parametric Bayes based on the structurally wrong model (Ricker)

```{r ricker-model}
```

Compute prior curves

```{r ricker-priors, dependson="common-priors"}
```

We define which parameters to keep track of, and set the initial values of
parameters in the transformed space used by the MCMC.

```{r ricker-mcmc, dependson="ricker-model"}
```

#### Convergence diagnostics for parametric bayes Ricker model

```{r ricker-traces, dependson="ricker-mcmc"}
```

```{r ricker-posteriors, dependson=c("ricker-traces", "ricker-priors")}
```

Reshape posteriors data, transform back, calculate mode and corresponding function.

```{r ricker-output, dependson=c("posterior-mode", "ricker-traces")}
```

## Myers Parametric Bayes

```{r myers-model}
```

```{r myers-priors}
```

```{r myers-mcmc-local, dependson="myers-model"}
jags.params=c("r0", "theta", "K", "stdQ")
jags.inits <- function(){
  list("r0"= 1.5 * rlnorm(1,0,.2),
       "K"=    8 * rlnorm(1,0,.1),
       "theta" = 2.5 * rlnorm(1,0,.1),
       "stdQ"= sqrt(0.2) * rlnorm(1,0,.1),
       .RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}
set.seed(12345)
myers_jags <- do.call(jags.parallel,
                      list(data=jags.data, inits=jags.inits,
    											 jags.params, n.chains=n.chains,
													 n.iter=n.iter, n.thin=n.thin,
                           n.burnin=n.burnin,
                           model.file="myers_process.bugs"))
recompile(myers_jags)
myers_jags <- do.call(autojags,
                      list(myers_jags, n.update=n.update,
                           n.iter=n.iter, n.thin = n.thin,
                           progress.bar="none"))
```

Convergence diagnostics for parametric bayes

```{r myers-traces, dependson="myers-mcmc-local"}
```

```{r myers-posteriors, dependson="myers-traces"}
```

```{r myers-output, dependson=c("posterior-mode", "myers-traces")}
```

### Phase-space diagram of the expected dynamics

```{r assemble-models, dependson=c("myers-output", "ricker-output", "allen-output", "gp-output", "mle-output")}
```

```{r Figure1, dependson="assemble-models"}
```

## Goodness of fit

This shows only the mean predictions.  For the Bayesian cases, we can instead loop over the posteriors of the parameters (or samples from the GP posterior) to get the distribution of such curves in each case.

```{r Figureb, dependson=c("assemble-models", "par-fns")}
```

## Optimal policies by value iteration

Compute the optimal policy under each model using stochastic dynamic programming. We begin with the policy based on the GP model,

```{r gp-opt, dependson="gp-output"}
```

Determine the optimal policy based on the allen and MLE models

```{r mle-opt, dependson="mle-output"}
```

Determine the optimal policy based on Bayesian Allen model

```{r allen-opt, dependson="allen-output"}
```

Bayesian Ricker

```{r ricker-opt, dependson="ricker-output"}
```

Bayesian Myers model

```{r myers-opt, dependson="myers-output"}
```

Assemble the data

```{r assemble-opt, dependson=c("gp-opt", "mle-opt", "allen-opt", "ricker-opt", "myers-opt")}
```

## Graph of the optimal policies

```{r Figure2, dependson="assemble-opt"}
```

## Simulate 100 realizations managed under each of the policies

```{r sims, dependson="assemble-opt"}
```

```{r Figure3, dependson="sims"}
```

```{r profits, dependson="sims"}
```

```{r totalprofits, dependson="profits"}
```

```{r deviances, dependson=c("posterior-mode", "myers-output", "ricker-output", "allen-output", "gp-output", "mle-output") }
```

repo_name: cboettig/nonparametric-bayes

repo_language: TeX

mime_type: text/plain


text:
# Comparison of Nonparametric Bayesian Gaussian Process estimates to standard the Parametric Bayesian approach

Plotting and knitr options, (can generally be ignored)

```{r plotting-options, cache=FALSE, message=FALSE, warning=FALSE, include=FALSE}
source("~/.knitr_defaults.R")
#opts_knit$set(upload.fun = socialR::flickr.url)
library(knitcitations)
library(nonparametricbayes)
opts_chunk$set(external=TRUE)
read_chunk("external-chunks.R")
```

```{r posterior-mode}
```

### Model and parameters

Uses the model derived in ` citet("10.1080/10236190412331335373")`, of a Ricker-like growth curve with an allee effect, defined in the pdgControl package,

```{r stateeq}
```

Various parameters defining noise dynamics, grid, and policy costs.

```{r sdp-pars, dependson="stateeq"}
```

### Sample Data

```{r obs, dependson="sdp-pars"}
```

## Maximum Likelihood

```{r mle, dependson="obs"}
```

Mean predictions

```{r mle-output, dependson="mle"}
```

## Non-parametric Bayes

```{r gp-priors}
```

Estimate the Gaussian Process (nonparametric Bayesian fit)

```{r gp, dependson=c("gp-priors", "obs")}
```

Show traces and posteriors against priors

```{r gp_traces_densities, dependson="gp"}
```

```{r gp-output, dependson="gp"}
```

Parametric Bayesian Models
---------------------------

We use the JAGS Gibbs sampler, a recent open source BUGS
implementation with an R interface that works on most platforms.
We initialize the usual MCMC parameters; see `?jags` for details.

All parametric Bayesian estimates use the following basic parameters for the JAGS MCMC:

```{r jags-setup}
```

We will use the same priors for process and observation noise in each model,

```{r common-priors}
```

### Parametric Bayes of correct (Allen) model

We initiate the MCMC chain (`init_p`) using the true values of the
parameters `p` from the simulation.  While impossible in real data, this
gives the parametric Bayesian approach the best chance at succeeding.
`y` is the timeseries (recall `obs` has the $x_t$, $x_{t+1}$ pairs)

The actual model is defined in a `model.file` that contains an R function
that is automatically translated into BUGS code by *R2WinBUGS*.  The file
defines the priors and the model. We write the file from R as follows:

```{r allen-model}
```

Write the priors into a list for later reference

```{r allen-priors, dependson="common-priors"}
```

We define which parameters to keep track of, and set the initial values of
parameters in the transformed space used by the MCMC.  We use logarithms
to maintain strictly positive values of parameters where appropriate.

```{r allen-mcmc, dependson=c("allen-model", "allen-pars", "jags-setup"), results="hide"}
```

#### Convergence diagnostics for Allen model

R notes: this strips classes from the `mcmc.list` object (so that we have list of matrices; objects that `reshape2::melt` can handle intelligently), and then combines chains into one array. In this array each parameter is given its value at each sample from the posterior (index) for each chain.

```{r allen-traces, dependson="allen-mcmc"}
```

```{r allen-posteriors, dependson=c("allen-traces", "allen-priors")}
```

Reshape the posterior parameter distribution data, transform back into original space, and calculate the mean parameters and mean function

```{r allen-output, dependson=c("posterior-mode", "allen-traces")}
```

## Parametric Bayes based on the structurally wrong model (Ricker)

```{r ricker-model}
```

Compute prior curves

```{r ricker-priors, dependson="common-priors"}
```

We define which parameters to keep track of, and set the initial values of
parameters in the transformed space used by the MCMC.

```{r ricker-mcmc, dependson="ricker-model"}
```

#### Convergence diagnostics for parametric bayes Ricker model

```{r ricker-traces, dependson="ricker-mcmc"}
```

```{r ricker-posteriors, dependson=c("ricker-traces", "ricker-priors")}
```

Reshape posteriors data, transform back, calculate mode and corresponding function.

```{r ricker-output, dependson=c("posterior-mode", "ricker-traces")}
```

## Myers Parametric Bayes

```{r myers-model}
```

```{r myers-priors}
```

```{r myers-mcmc, dependson="myers-model"}
```

Convergence diagnostics for parametric bayes

```{r myers-traces, dependson="myers-mcmc"}
```

```{r myers-posteriors, dependson="myers-traces"}
```

```{r myers-output, dependson=c("posterior-mode", "myers-traces")}
```

### Phase-space diagram of the expected dynamics

```{r assemble-models, dependson=c("myers-output", "ricker-output", "allen-output", "gp-output", "mle-output")}
```

```{r Figure1, dependson="assemble-models"}
```

## Goodness of fit

This shows only the mean predictions.  For the Bayesian cases, we can instead loop over the posteriors of the parameters (or samples from the GP posterior) to get the distribution of such curves in each case.

```{r Figureb, dependson=c("assemble-models", "par-fns")}
```

## Optimal policies by value iteration

Compute the optimal policy under each model using stochastic dynamic programming. We begin with the policy based on the GP model,

```{r gp-opt, dependson="gp-output"}
```

Determine the optimal policy based on the allen and MLE models

```{r mle-opt, dependson="mle-output"}
```

Determine the optimal policy based on Bayesian Allen model

```{r allen-opt, dependson="allen-output"}
```

Bayesian Ricker

```{r ricker-opt, dependson="ricker-output"}
```

Bayesian Myers model

```{r myers-opt, dependson="myers-output"}
```

Assemble the data

```{r assemble-opt, dependson=c("gp-opt", "mle-opt", "allen-opt", "ricker-opt", "myers-opt")}
```

## Graph of the optimal policies

```{r Figure2, dependson="assemble-opt"}
```

## Simulate 100 realizations managed under each of the policies

```{r sims, dependson="assemble-opt"}
```

```{r Figure3, dependson="sims"}
```

```{r profits, dependson="sims"}
```

```{r totalprofits, dependson="profits"}
```

```{r deviances, dependson=c("posterior-mode", "myers-output", "ricker-output", "allen-output", "gp-output", "mle-output") }
```


text:
#' Plot a Gaussian Process
#'
#' Plot the mean (line) and 2 standard deviations (ribbon)
#' @param x a gpfit object from \code{\link{gp_fit}}
#' @param true an optional data.frame with columns "x" and "y" containing target curve.
#' @param ... additional arguments (not implemented)
#' @return a ggplot object of the plot
#' @export
plot.gpfit <- function(x, true, ...){
  dat <- data.frame(x=x$X, y=(x$Ef), ymin=(x$Ef-2*sqrt(diag(x$Cf))), ymax=(x$Ef+2*sqrt(diag(x$Cf))))
  p <- ggplot(dat) +
    geom_ribbon(aes(x=x,y=y, ymin=ymin, ymax=ymax), fill="grey80") + # Var
    geom_line(aes(x=x,y=y), size=1) + #MEAN
    geom_point(data=gp$obs,aes(x=x,y=y)) +
    labs(title=paste("llik =", prettyNum(gp$llik)))
  if(!is.null(true))
    p <- p + geom_line(data = true, aes(x, y), col = "red", lty = 2)
  p
}


text:
#' Determine the transition matrix at each h given a Guassian process
#'
#' @param Ef GP mean
#' @param V GP vector of variances
#' @param x_grid x grid over which to define transition matrices
#' @param h_grid the harvest grid over which to compute.  Assumed to match the stock grid x_grid if not given.
#' @return a list of matrices corresponding to the state transition matrix at each harvest level h
#' @export
gp_transition_matrix <- function(Ef, V, x_grid, h_grid=NULL){
  if(is.null(h_grid))
    h_grid <- x_grid
  if(is.matrix(Ef)){ # loop over posteriors
    n_samples <- dim(Ef)[2]
    loop_over_h <-function(h){
      # Average across matrices in list, and then normalize
      F_s <- gp_calc_F(h, Ef[,1], V[,1], x_grid)
      for(i in 2:n_samples){
        F_s <- F_s + gp_calc_F(h, Ef[,i], V[,i], x_grid)
      }
      F_s <- rownorm(F_s)
    }
    out <- lapply(h_grid, loop_over_h)
  } else {
    out <- lapply(h_grid, gp_calc_F, Ef, V, x_grid)
  }
  out
}

gp_calc_F <- function(h, Ef, V, x_grid){
  j <- which.min(abs(x_grid - h))
  indices <- 1:length(Ef)
  shift <- pmax(indices-j,1)
  mu <- Ef[shift]
  mu[shift == 1] <- 0
  F_ <- t(sapply(1:length(mu), function(i){
    if(snap_to_grid(mu[i],x_grid) < x_grid[2]){ #
      out <- numeric(length(x_grid))
      out[1] <- 1
      out
    } else {
      out <- dnorm(x_grid, mu[i], sqrt(V[i]))
    }
  }))
  F_ <- rownorm(F_)
}

#' calculate transition matrix at each harvest level based on parameterized function
#'
#' @param f the transtion function, a function of state x, harvest h, and parameters p
#' @param p the parameters p of the transition function
#' @param x_grid the grid over which to solve
#' @param h_grid the harvest grid over which to compute.  Assumed to match the stock grid gp$X if not given.
#' @param sigma_g the scale of the noise
#' @return a list of matrices corresponding to the state transition matrix at each harvest level h
#' @export
f_transition_matrix <- function(f, p, x_grid, h_grid = NULL, sigma_g){
  if(is.null(h_grid))
    h_grid <- x_grid
  lapply(h_grid, par_F, f, p, x_grid, sigma_g)
}

#' calculate transition matrix based on parameterized function
#'
#' @param h the harvest over which to compute.
#' @param f the transtion function, a function of state x, harvest h, and parameters p
#' @param p the parameters p of the transition function
#' @param x_grid the grid over which to solve
#' @param sigma_g the scale of the noise
#' @return a list of matrices corresponding to the state transition matrix at each harvest level h
#' @export
par_F <- function(h, f, p, x_grid, sigma_g){
  mu <- sapply(x_grid, f, h, p)
  F_true <- t(sapply(mu, function(m){
    if(snap_to_grid(m,x_grid) < x_grid[2]){ #
      out <- numeric(length(x_grid))
      out[1] <- 1
      out
    } else {
      out <- dlnorm(x_grid/m, 0, sigma_g)
    }
  }))
  F_true <- rownorm(F_true)
}

snap_to_grid <- function(x, grid) sapply(x, function(x) grid[which.min(abs(grid - x))])

# internal helper function
rownorm <- function(M)
  t(apply(M, 1, function(x){
    if(sum(x)>0){
      x/sum(x)
    } else {
      out <- numeric(length(x))
      out[1] <- 1
      out
    }
  }))


text:
#' Simulate observed data for training the GP
#'
#' Simulate observed data for training the GP
#' @param Xo the initial population size
#' @param z_g a function of no arguments providing the random multiplicative growth noise
#' @param f a function of (x, h, p) providing the growth in population x under harvest h given pars p
#' @param p the parameters of function f
#' @param Tobs the number of timesteps on which to observe the system
#' @param seed the random number seed for the simulation
#' @param nz the number of observations at 0 (a hack, should just condiition on 0,0)
#' @param harvest a sequence of harvest values that occur durning the simulation
#' @param plotdata logical, plot data?
#' @details varying harvest values allow the system to explore the state space, making for better training data.
#' @noRd
sim_obs <- function(Xo, z_g, f, p, Tobs = 35, seed = 1, nz = 10,
                    harvest = sort(rep(seq(0, .5, length=7), 5)),
                    plotdata=FALSE){
  x <- numeric(Tobs)
  x[1] <- Xo
  set.seed(seed)
  for(t in 1:(Tobs-1))
    x[t+1] = z_g() * f(x[t], h=harvest[t], p=p)
  if(plotdata)
    plot(x)
  obs <- data.frame(x = c(rep(0,nz),
                          pmax(rep(0,Tobs-1), x[1:(Tobs-1)] -
                                 harvest[1:Tobs-1])),
                    y = c(rep(0,nz),
                          x[2:Tobs]))
  #plot(obs$x, obs$y)
  obs
}

#' MLE estimate of parameters given by function f
#'
#' MLE estimate of parameters given by function f
#' @param obs the observed data, two columns giving x_t, x_t+1 respectively
#' @return a list with the f given, the MLE estimated parameters, and estimated noise level
#' @noRd
par_est_allee <- function(obs, f, p,
                          init = c(p[1]+1, p[2]-1, p[3]-2,
                                   s = sigma_g + abs(rnorm(1,0,.1)))
                          ){
  estf <- function(p){
    mu <- f(obs$x,0,p)
    -sum(dlnorm(obs$y, log(mu), p[4]), log=TRUE)
  }
  par = init
  o <- optim(par, estf, method="L", lower=c(1e-3,1e-3,1e-3, 1e-3))
  f_alt <- f
  p_alt <- c(as.numeric(o$par[1]), as.numeric(o$par[2]), as.numeric(o$par[3]))
  sigma_g_alt <- as.numeric(o$par[4])
  list(f = f_alt, p = p_alt, sigma_g = sigma_g_alt)
}

#' MLE estimate of parameters given by function f
#' @param obs the observed data, two columns giving x_t, x_t+1 respectively
#' @return a list with the f given, the MLE estimated parameters, and estimated noise level
#' @noRd
par_est <- function(obs,  init = c(r=1.5, K=mean(obs$x), s=1)){
  estf <- function(p){
    mu <- log(obs$x) + p["r"]*(1-obs$x/p["K"])
    llik <- -sum(dlnorm(obs$y, mu, p["s"]), log=TRUE)
    if(!is.numeric(llik) | is.nan(llik) | !(llik < Inf)){
      warning("possible error in llik")
      llik <- 1e30
    }
    llik
  }
  o <- optim(par = init, estf, method="L", lower=c(1e-3,1e-3,1e-3))
  f_alt <- Ricker
  p_alt <- c(o$par['r'], o$par['K'])
  sigma_g_alt <- o$par['s']
  list(f_alt = f_alt, p_alt = p_alt, sigma_g_alt = sigma_g_alt)
}

#' helper function to determine the optimal policies of each model

#' helper function to determine the optimal policies of each model
#' @import pdgControl
#' @noRd
optimal_policy <- function(gp, f, f_est, f_alt, p, p_est, p_alt, x_grid, h_grid, sigma_g, sigma_g_est, sigma_g_alt, delta, xT, profit, reward, OptTime){
  MaxT = 1000
  matrices_gp <- gp_transition_matrix(gp$ZZ.km, gp$ZZ.ks2, x_grid, h_grid)
  opt_gp <- value_iteration(matrices_gp, x_grid, h_grid, MaxT, xT, profit, delta, reward)
  matrices_true <- f_transition_matrix(f, p, x_grid, h_grid, sigma_g)
  opt_true <- value_iteration(matrices_true, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
  matrices_estimated <- f_transition_matrix(f_est, p_est, x_grid, h_grid, sigma_g_est)
  opt_estimated <- value_iteration(matrices_estimated, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)
  matrices_alt <- f_transition_matrix(f_alt, p_alt, x_grid, h_grid, sigma_g_alt)
  opt_alt <- value_iteration(matrices_alt, x_grid, h_grid, OptTime=MaxT, xT, profit, delta=delta)

  list(gp_D = opt_gp$D, true_D = opt_true$D, est_D = opt_estimated$D, alt_D = opt_alt$D)
}

#' helper function to simulate the optimal policies of each model
#' @import pdgControl data.table reshape2
#' @noRd
simulate_opt <- function(OPT, f, p, x_grid, h_grid, x0, z_g, profit, OptTime){
  set.seed(1)
  sim_gp <- lapply(1:100, function(i)
    ForwardSimulate(f, p, x_grid, h_grid, x0, OPT$gp_D, z_g, profit=profit, OptTime=OptTime))
  set.seed(1)
  sim_true <- lapply(1:100, function(i)
    ForwardSimulate(f, p, x_grid, h_grid, x0, OPT$true_D, z_g, profit=profit, OptTime=OptTime))
  set.seed(1)
  sim_est <- lapply(1:100, function(i)
    ForwardSimulate(f, p, x_grid, h_grid, x0, OPT$est_D, z_g, profit=profit, OptTime=OptTime))
  set.seed(1)
  sim_alt <- lapply(1:100, function(i)
    ForwardSimulate(f, p, x_grid, h_grid, x0, OPT$alt_D, z_g, profit=profit, OptTime=OptTime))

  dat <- list(GP = sim_gp, Parametric = sim_est,
              True = sim_true, Structural = sim_alt)
  dat <- melt(dat, id=names(dat[[1]][[1]]))
  dt <- data.table(dat)
  setnames(dt, c("L1", "L2"), c("method", "reps"))
  setkey(dt, method) # chance the ordering
  dt[c("GP", "Parametric", "True", "Structural")]
}

#' helper function to draw the plots of fish stock and harvest dynamics
#' @param dt the data.table from simulate_opt
#' @import ggplot2
#' @return plots of stock and harvest dynamics over time and replicates
#' @noRd
sim_plots <- function(dt, seed=1, harvest=FALSE){
  fish_plot <- ggplot(dt) +
    geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1) +
    scale_colour_manual(values=cbPalette, guide = guide_legend(override.aes = list(alpha = 1))) +
    scale_colour_manual(values=c("GP" = cbPalette[1], "Parametric" = cbPalette[2], "True" = cbPalette[3], "Structural" = cbPalette[4]),
                        guide = guide_legend(override.aes = list(alpha = 1))) +
    ggtitle(paste("seed", seed))
  print(fish_plot)
  harvest_plot <- ggplot(dt) +
    geom_line(aes(time, harvest, group=interaction(reps,method), color=method), alpha=.1) +
    scale_colour_manual(values=c("GP" = cbPalette[1], "Parametric" = cbPalette[2], "True" = cbPalette[3], "Structural" = cbPalette[4]),
                        guide = guide_legend(override.aes = list(alpha = 1))) +
    ggtitle(paste("seed", seed))
  if(harvest)
    print(harvest_plot)
}

#' function to provide summary stats or the profits
#' @param dt the data table from simulate_opt
#' @return a table of mean and variance in profit over the replicates, by method
#' @import data.table
#' @noRd
profits_stats <- function(dt){
  profits <- dt[, sum(profit), by = c("reps", "method")]
  means <- profits[, mean(V1), by = method]
  sds <- profits[, sd(V1), by = method]
  yield <- cbind(means, sd = sds$V1)
  yield
}

#' plot the gaussian process, true model, and fitted parametric model(s)
#' @import ggplot2 reshape2
#' @noRd
gp_plot <- function(gp, f, p, f_est, p_est, f_alt, p_alt, x_grid, obs, seed){
  tgp_dat <-
    data.frame(  x = gp$XX[[1]],
                 y = gp$ZZ.km,
                 ymin = gp$ZZ.km - 1.96 * sqrt(gp$ZZ.ks2),
                 ymax = gp$ZZ.km + 1.96 * sqrt(gp$ZZ.ks2),
                 ymin2 = gp$ZZ.mean - 1.96 * sqrt(gp$ZZ.vark),
                 ymax2 = gp$ZZ.mean + 1.96 * sqrt(gp$ZZ.vark))
  true <- sapply(x_grid, f, 0, p)
  alt <- sapply(x_grid, f_alt, 0, p_alt)
  est <- sapply(x_grid, f_est, 0, p_est)
  models <- data.frame(x=x_grid, GP=tgp_dat$y, Parametric=est, True=true, Structural=alt)
  models <- melt(models, id="x")
  names(models) <- c("x", "method", "value")
  plot_gp <- ggplot(tgp_dat) + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") +
    geom_ribbon(aes(x,y,ymin=ymin2,ymax=ymax2), fill="gray60") +
    geom_line(data=models, aes(x, value, col=method), lwd=2, alpha=0.8) +
    geom_point(data=obs, aes(x,y), alpha=0.8) +
    xlab(expression(X[t])) + ylab(expression(X[t+1])) +
    scale_colour_manual(values=cbPalette) +
    ggtitle(paste("seed", seed))
  print(plot_gp)
}

#' plot the gp posteriors
#' @param gp a gp fit
#' @param priors prior distribution functions so that they can be shown on plot as well
#' @return the posteriors plot as histogram, with priors overlaid
#' @noRd
posteriors_plot <- function(gp, priors){
  hyperparameters <- c("index", "s2", "tau2", "beta0", "nug", "d", "ldetK")
  posteriors <- melt(gp$trace$XX[[1]][,hyperparameters], id="index")
  prior_curves <- ddply(posteriors, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = priors[[dd$variable[1]]](grid))
  })
  plot_posteriors <- ggplot(posteriors) +
    #geom_density(aes(value), lwd=2) +
    geom_histogram(aes(x=value, y=..density..), alpha=0.7) +
    geom_line(data=prior_curves, aes(x=value, y=density), col="red", lwd=2) +
    facet_wrap(~ variable, scale="free")
  print(plot_posteriors)
}

#' policy plot
#' @param x_grid stock grid
#' @param gpD gaussian process policy
#' @param estD an estimated policy (e.g. parametric uncertainty only)
#' @param trueD the true optimal policy
#' @param altD an alternative (structurally incorrect) policy
#' @return the policy plot
#' @noRd
plot_policies <- function(x_grid, gpD, estD, trueD, altD){
  policies <- melt(data.frame(stock=x_grid,
                              GP = x_grid[gpD],
                              Parametric = x_grid[estD],
                              True = x_grid[trueD],
                              Structural = x_grid[altD]),
                   id="stock")
  names(policies) <- c("stock", "method", "value")
  policy_plot <- ggplot(policies, aes(stock, stock - value, color=method)) +
    geom_line(lwd=2, alpha=0.8) +
    xlab("stock size") + ylab("escapement")  +
    scale_colour_manual(values=cbPalette)
  print(policy_plot)
}


text: cbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")


text:
# edit this file for a demo

#' gp_fit fits a Gaussian process to the observed data
#'
#' API not yet stable! Arguments to this function will probably be updated to be much more generic
#' @param obs a data frame of observations with columns obs$x and obs$y
#' @param X the desired points over which to predict
#' @param pars a named numeric specifying "sigma_n" for the (additive) noise and "l" for the covariance length scale
#' @param method select the method to use
#' @param out_var optionally force kernlab method to generate Cf using a different variance than sigma_n^2 assumed in the fit
#' @param fit logical, argument for the kernlab::gausspr method only
#' @param ... potential additional arguments, unimplemented
#' @return a list with items "mu", the expected Y values at X (mean of the posterior Gaussian process),
#' Sigma, the covariance matrix for the posterior Gaussian process, and
#' "loglik", the log likelihood of observering the given data under the process, marginalized over the prior
#' @details so far treates the prior as mean 0 and covariance given by cov(X)
#' @import kernlab
#' @export
# @examples
# X <- seq(-5,5,len=50)
# obs <- data.frame(x = c(-4, -3, -1,  0,  2),
#                   y = c(-2,  0,  1,  2, -1))
# gp <- gp_fit(obs, X, c(sigma_n=0.8, l=1))
#
# df <- data.frame(x=X, y=gp$mu, ymin=(gp$mu-2*sqrt((gp$var))), ymax=(gp$mu+2*sqrt((gp$var))))
# require(ggplot2)
# ggplot(df)  + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") + geom_line(aes(x,y)) + geom_point(data=obs, aes(x,y))

gp_fit <- function(obs, X, pars=c(sigma_n=1, tau=1, l=1), method=c("direct", "sequential", "kernlab", "cholesky", "conditional"), fit = TRUE, out_var=NULL, ...){

  method <- match.arg(method)
  sigma_n <- pars["sigma_n"]
  l <- pars["l"]
  tau <- pars["tau"]
  llik <- NA

  ## Parameterization-specific
  SE <- function(Xi,Xj, l=l) tau * exp(-0.5 * (Xi - Xj) ^ 2 / l ^ 2)
  cov <- function(X, Y) outer(X, Y, SE, l)
  n <- length(obs$x)
  K <- cov(obs$x, obs$x)
  I <-  diag(1, n)

  ## Cholesky simultaneous method.
  if(method=="cholesky"){
    L <- t(chol(K + sigma_n ^ 2 * I))
    alpha <- solve(t(L), solve(L, obs$y))
    loglik <- -.5 * t(obs$y) %*% alpha - sum(log(diag(L))) - n * log(2 * pi) / 2
    k_star <- cov(obs$x, X)
    mu <- t(k_star) %*% alpha
    v <- solve(L, k_star)
    var <- diag( cov(X,X) - t(v) %*% v )
  }

  ## Direct method
  else if(method=="direct"){
    cov_xx_inv <- solve(K + sigma_n ^ 2 * I)
    Ef <- cov(X, obs$x) %*% cov_xx_inv %*% obs$y
    Cf <- cov(X, X) - cov(X, obs$x)  %*% cov_xx_inv %*% cov(obs$x, X)
    llik <- -.5 * t(obs$y) %*% cov_xx_inv %*% obs$y - 0.5 * log(det(cov_xx_inv)) - n * log(2 * pi) / 2
  }

  ## Direct method, conditional on passing through 0,0 (quick hack)
  else if(method=="conditional"){
    obs <- data.frame(x=c(0, obs$x), y = c(0, obs$y))
    n <- length(obs$x)
    K <- cov(obs$x, obs$x)
    sigmaI <-  sigma_n ^ 2 * diag(1, n)
    sigmaI[1,1] <- 1e-5 # no noise about zero.
    cov_xx_inv <- solve(K + sigmaI)
    Ef <- cov(X, obs$x) %*% cov_xx_inv %*% obs$y
    Cf <- cov(X, X) - cov(X, obs$x)  %*% cov_xx_inv %*% cov(obs$x, X)
    llik <- -.5 * t(obs$y) %*% cov_xx_inv %*% obs$y - 0.5 * log(det(cov_xx_inv)) - n * log(2 * pi) / 2
  }

  ### Sequential method, using function recursion (terribly slow for more than a few data points)
  else if(method=="sequential"){
    mmult <- function(x,y){
      if(length(x) == 1){
        x <- as.numeric(x)
        return(x * y)
      }
      if(length(y) == 1){
        y <- as.numeric(y)
        return(x * y)
      }
      return(x %*% y)
    }

    mu <- numeric(length(obs$y))
    y <- obs$y
    x <- obs$x
    C_seq <- function(X, X_prime, i){
      if(i <= 1)
        cov(X, X_prime) - mmult(cov(X,x[i]), cov(x[i], X_prime)) / as.numeric( cov(x[i], x[i]) + sigma_n^2)
      else
        C_seq(X, X_prime,   i-1) - mmult(C_seq(X,x[i],   i-1), C_seq(x[i], X_prime,   i-1)) / as.numeric( C_seq(x[i], x[i],   i-1)  + sigma_n^2  )
    }
    mu_seq <- function(X, i){
      if(i <= 1)
        cov(x[i], X) * (y[i]-mu[i]) / as.numeric( cov(x[i], x[i]) + sigma_n^2)
      else
        mu_seq(X, i-1) + C_seq(x[i], X,  i-1) * (y[i]-mu[i]) / as.numeric( C_seq(x[i], x[i], i-1)  + sigma_n^2 )
    }
    Ef <- t(mu_seq(X, length(obs$x)))
    Cf <- C_seq(X, X, length(obs$x))
  }

  ## Kernlab method
  else if(method=="kernlab"){
  var <- sigma_n^2
  if(!fit)
    gp <- gausspr(obs$x, obs$y, kernel="rbfdot", kpar=list(sigma=1/(2*l^2)), fit=FALSE, scaled=FALSE, var=var) # fixed length-scale
  else
    gp <- gausspr(obs$x, obs$y, kernel="rbfdot", scaled=FALSE, var = var)
  Ef <- predict(gp, X)

  if(!is.null(out_var)) # permit an external variance
    var <- out_var

  ##Manually get the covariance, $K(x_*, x_*) - K(x_*, x)(K(x, x) - \sigma_n^2\mathbb{I})^{-1}K(x,x_*)$.
  Inv <- solve(kernelMatrix(kernelf(gp), xmatrix(gp)) + diag(rep(var, length = length(xmatrix(gp)))))
  Cf <- kernelMatrix(kernelf(gp), as.matrix(X)) -
    kernelMult(kernelf(gp), as.matrix(X), xmatrix(gp),  Inv) %*%
    kernelMatrix(kernelf(gp), xmatrix(gp), as.matrix(X))
  llik <- -.5 * t(obs$y) %*% Inv %*% obs$y - 0.5 * log(det(Inv)) - n * log(2 * pi) / 2

  }

  out <- list(Ef = Ef, Cf = Cf, llik = llik, obs = obs, X = X, K = K, pars = pars)
  class(out) = "gpfit"
  out
  }

update.gpfit <- function(gp, obs, sigma){

}


text:
#' Basic regression in Gaussian processes
#'
#'
#' @param x Observed x values, (vector or matrix with columns for each dimension of data)
#' @param y Vector of observed y values in the training data
#' @param init_pars the initial guesses for lengthscale l and process noise sigma_n
#' @param n iterations of the metropolis algorithm
#' @param d.p parameters for the length-scale prior, as an inverse Gamma distribution
#' @param s2.p parameters for the noise prior, as an inverse Gamma distribution
#' @details Currently assumes the covariance function.  By default we will use
#' the squared exponential (also called radial basis or Gaussian,
#' though it is not this that gives Gaussian process it's name;
#' any covariance function would do) as the covariance function,
#' \deqn{\operatorname{cov}(X_i, X_j) = \exp\left(-\frac{(X_i - X_j)^2}{2 \ell ^ 2}\right)}{}
#' @return the MCMC output, as constructed by \link{metrop} from the mcmc package.
#' @export
#' @import mcmc
gp_mcmc <- function(x, y, init_pars = c(l=1, sigma.n=1), n = 1e4, d.p = c(5,5), s2.p = c(5,5)){

  lpriors <- function(pars){
    prior <- unname(
      dgamma(exp(pars[1]), d.p[1], scale = d.p[2]) *
      dgamma(exp(pars[2]), s2.p[1], s2.p[2])
    )
    log(prior)
  }
  posterior <- function(pars, x, y){
    l <- exp(pars[1])
    sigma.n <- exp(pars[2])

    cov <- function(X, Y) outer(X, Y, SE, l)
    I <- diag(1, length(x))
    K <- cov(x, x)

    loglik <- - 0.5 * t(y) %*% solve(K + sigma.n^2 * I) %*% y -
      log(det(K + sigma.n^2*I)) -
      length(y) * log(2 * pi) / 2

    loglik + lpriors(pars)
  }

  out <- metrop(posterior, log(init_pars), n, x = x, y = y)

  # assemble output
  obs = list(x=x,y=y)
  out$d.p <- d.p
  out$s2.p <- s2.p
  out$obs <- obs
  out
}

#' predict the expected values and posterior distributions of the Gaussian Process
#'
#' @param gp a fit of the gaussian process from gp_mcmc
#' @param x_predict the values at which we desire predictions
#' @param burnin length of sequence to discard as a transient
#' @param thin frequency of sub-sampling (make posterior distribution smaller if necessary)
#' @param covs if TRUE, return covariances instead of variances (set high thinning as this is memory intensive)
#' @import reshape2
#' @export
gp_predict <- function(gp, x_predict, burnin=0, thin=1, covs=FALSE){

  postdist <- cbind(index=1:gp$nbatch, as.data.frame(exp(gp$batch)))
  s <- seq(burnin+1, gp$nbatch, by=thin)
  postdist <- postdist[s,]
  names(postdist) <- c("index", names(gp$initial))
  indices <- 1:dim(postdist)[1]
  obs <- gp$obs

  out <- lapply(indices, function(i){
    l <- postdist[i, "l"]
    sigma.n <- postdist[i, "sigma.n"]

    cov <- function(X, Y) outer(X, Y, SE, l)
    cov_xx_inv <- solve(cov(obs$x, obs$x) + sigma.n^2 * diag(1, length(obs$x)))
    Ef <- cov(x_predict, obs$x) %*% cov_xx_inv %*% obs$y
    Cf <- cov(x_predict, x_predict) - cov(x_predict, obs$x)  %*% cov_xx_inv %*% cov(obs$x, x_predict)
    list(Ef = Ef, Cf = Cf, Vf = diag(Cf))
  })

  Ef_posterior <- sapply(out, `[[`, "Ef")
  Cf_posterior <- sapply(out, `[[`, "Cf")
  Vf_posterior <- sapply(out, `[[`, "Vf")

  E_Ef <- rowMeans(Ef_posterior)
  E_Cf <- matrix( apply(Cf_posterior, 1, sum) / dim(Cf_posterior)[2], ncol = sqrt(dim(Cf_posterior)[1]) )
  E_Vf <- diag(E_Cf) # same as rowMeans(Vf_posterior)

  # list format better for return
  Cf_posterior <- lapply(out, `[[`, "Cf")

  list(Ef_posterior = Ef_posterior,
       Vf_posterior = Vf_posterior,
       Cf_posterior = Cf_posterior,
       E_Ef = E_Ef, E_Cf = E_Cf, E_Vf = E_Vf)

}

#' Summary plots showing the trace and posteriors for the gp_mcmc estimates
#'
#' @param gp a fit of the gaussian process from gp_mcmc
#' @param burnin length of sequence to discard as a transient
#' @param thin frequency of sub-sampling (make posterior distribution smaller if necessary)
#' @return two ggplot2 objects, one plotting the trace and one
#' plotting the posteriors in black with priors overlaid in red.
#' @import reshape2
#' @import ggplot2
#' @import plyr
#' @export
summary_gp_mcmc <- function(gp, burnin=0, thin=1){

  postdist <- cbind(index=1:gp$nbatch, as.data.frame(exp(gp$batch)))
  s <- seq(burnin+1, gp$nbatch, by=thin)
  postdist <- postdist[s,]
  names(postdist) <- c("index", names(gp$initial))

  # TRACES
  df <- melt(postdist, id="index")
  traces_plot <-
    ggplot(df) + geom_line(aes(index, value)) +
    facet_wrap(~ variable, scales="free", ncol=1)

  s2.p <- gp$s2.p
  d.p <- gp$d.p
  s2_prior <- function(x) dgamma(x, s2.p[1], s2.p[2])
  d_prior <- function(x) dgamma(x, d.p[1], scale = d.p[2])
  prior_fns <- list(l = d_prior, sigma.n = s2_prior)

  prior_curves <- ddply(df, "variable", function(dd){
    grid <- seq(min(dd$value), max(dd$value), length = 100)
    data.frame(value = grid, density = prior_fns[[dd$variable[1]]](grid))
  })

  # Posteriors (easier to read than histograms)
  posteriors_plot <- ggplot(df, aes(value)) +
    stat_density(geom="path", position="identity", alpha=0.7) +
    geom_line(data=prior_curves, aes(x=value, y=density), col="red") +
    facet_wrap(~ variable, scales="free", ncol=2)

 # print(traces_plot)
#  print(posteriors_plot)

  out <- list(traces_plot = traces_plot, posteriors_plot = posteriors_plot)
  invisible(out)
}

# Helper function: The default covariance function
SE <- function(Xi,Xj, l) exp(-0.5 * (Xi - Xj) ^ 2 / l ^ 2)

# In general we aren't interested in drawing from the prior, but want to include information from the data as well.  We want the joint distribution of the observed values and the prior is:
# $$\begin{pmatrix} y_{\textrm{obs}} \\ y_{\textrm{pred}} \end{pmatrix} \sim \mathcal{N}\left( \mathbf{0}, \begin{bmatrix} cov(X_o,X_o) & cov(X_o, X_p) \\ cov(X_p,X_o) & cov(X_p, X_p) \end{bmatrix} \right)$$
## Additive noise
## In general the model may have process error, and rather than
## observe the deterministic map $f(x)$ we only observe $y = f(x) + \varepsilon$.
## Let us assume for the moment that $\varepsilon$ are independent, normally
## distributed random variables with variance $\sigma_n^2$.
##
## The likelihood is given by
# $$\log(p(y | X)) = -\tfrac{1}{2} \mathbf{y}^T (K + \sigma_n^2 \mathbf{I})^{-1} y  - \tfrac{1}{2} \log\left| K + \sigma_n^2 \mathbf{I} \right| - \tfrac{n}{2}\log 2 \pi$$


text:
#!/bin/sh
exec "${REALGCC:-gcc}" "$@" -specs "musl-gcc.specs"

repo_language: C

mime_type: text/x-c


text:
This is [Lua](http://www.lua.org/) + [SQLite](http://sqlite.org/) + [musl](http://www.musl-libc.org/) running on x86.

An x86 CPU capable of running in 64-bit mode and a graphics chipset supporting the VESA
BIOS Extensions are required.

Use [luakernel.iso](bin/luakernel.iso) with your favorite x86 system.

![](screenshot.png?raw=true)


text:
OUTPUT_FORMAT(elf64-x86-64)
ENTRY(start)
SECTIONS
{
  . = 0x100000;
  .text :
  {
    text_start = .;
    *(.text)
    text_end = .;
  }
  .data :
  {
    data_start = .;
    *(.data)
    data_end = .;
  }
  .bss :
  {
    bss_start = .;
    *(.bss)
    bss_end = .;
  }
  .heap :
  {
    . = ALIGN(16);
    heap_start = .;
  }
}


text:
%rename cpp_options old_cpp_options

*cpp_options:
-nostdinc -isystem dep/musl/include -isystem include%s %(old_cpp_options)

*cc1:
%(cc1_cpu) -nostdinc -isystem dep/musl/include -isystem include%s

*link_libgcc:
-L dep/musl/lib/

*libgcc:
libgcc.a%s %:if-exists(libgcc_eh.a%s)

*startfile:
%{!shared: dep/musl/lib/%{pie:S}crt1.o} dep/musl/lib/crti.o %{shared|pie:crtbeginS.o%s;:crtbegin.o%s}

*endfile:
%{shared|pie:crtendS.o%s;:crtend.o%s} dep/musl/lib/crtn.o

*link:
-dynamic-linker /lib/ld-musl-x86_64.so.1 -nostdlib %{shared:-shared} %{static:-static} %{rdynamic:-export-dynamic}

*esp_link:

*esp_options:

*esp_cpp_options:


text:
CC0 1.0 Universal

Statement of Purpose

The laws of most jurisdictions throughout the world automatically confer
exclusive Copyright and Related Rights (defined below) upon the creator and
subsequent owner(s) (each and all, an "owner") of an original work of
authorship and/or a database (each, a "Work").

Certain owners wish to permanently relinquish those rights to a Work for the
purpose of contributing to a commons of creative, cultural and scientific
works ("Commons") that the public can reliably and without fear of later
claims of infringement build upon, modify, incorporate in other works, reuse
and redistribute as freely as possible in any form whatsoever and for any
purposes, including without limitation commercial purposes. These owners may
contribute to the Commons to promote the ideal of a free culture and the
further production of creative, cultural and scientific works, or to gain
reputation or greater distribution for their Work in part through the use and
efforts of others.

For these and/or other purposes and motivations, and without any expectation
of additional consideration or compensation, the person associating CC0 with a
Work (the "Affirmer"), to the extent that he or she is an owner of Copyright
and Related Rights in the Work, voluntarily elects to apply CC0 to the Work
and publicly distribute the Work under its terms, with knowledge of his or her
Copyright and Related Rights in the Work and the meaning and intended legal
effect of CC0 on those rights.

1. Copyright and Related Rights. A Work made available under CC0 may be
protected by copyright and related or neighboring rights ("Copyright and
Related Rights"). Copyright and Related Rights include, but are not limited
to, the following:

  i. the right to reproduce, adapt, distribute, perform, display, communicate,
  and translate a Work;

  ii. moral rights retained by the original author(s) and/or performer(s);

  iii. publicity and privacy rights pertaining to a person's image or likeness
  depicted in a Work;

  iv. rights protecting against unfair competition in regards to a Work,
  subject to the limitations in paragraph 4(a), below;

  v. rights protecting the extraction, dissemination, use and reuse of data in
  a Work;

  vi. database rights (such as those arising under Directive 96/9/EC of the
  European Parliament and of the Council of 11 March 1996 on the legal
  protection of databases, and under any national implementation thereof,
  including any amended or successor version of such directive); and

  vii. other similar, equivalent or corresponding rights throughout the world
  based on applicable law or treaty, and any national implementations thereof.

2. Waiver. To the greatest extent permitted by, but not in contravention of,
applicable law, Affirmer hereby overtly, fully, permanently, irrevocably and
unconditionally waives, abandons, and surrenders all of Affirmer's Copyright
and Related Rights and associated claims and causes of action, whether now
known or unknown (including existing as well as future claims and causes of
action), in the Work (i) in all territories worldwide, (ii) for the maximum
duration provided by applicable law or treaty (including future time
extensions), (iii) in any current or future medium and for any number of
copies, and (iv) for any purpose whatsoever, including without limitation
commercial, advertising or promotional purposes (the "Waiver"). Affirmer makes
the Waiver for the benefit of each member of the public at large and to the
detriment of Affirmer's heirs and successors, fully intending that such Waiver
shall not be subject to revocation, rescission, cancellation, termination, or
any other legal or equitable action to disrupt the quiet enjoyment of the Work
by the public as contemplated by Affirmer's express Statement of Purpose.

3. Public License Fallback. Should any part of the Waiver for any reason be
judged legally invalid or ineffective under applicable law, then the Waiver
shall be preserved to the maximum extent permitted taking into account
Affirmer's express Statement of Purpose. In addition, to the extent the Waiver
is so judged Affirmer hereby grants to each affected person a royalty-free,
non transferable, non sublicensable, non exclusive, irrevocable and
unconditional license to exercise Affirmer's Copyright and Related Rights in
the Work (i) in all territories worldwide, (ii) for the maximum duration
provided by applicable law or treaty (including future time extensions), (iii)
in any current or future medium and for any number of copies, and (iv) for any
purpose whatsoever, including without limitation commercial, advertising or
promotional purposes (the "License"). The License shall be deemed effective as
of the date CC0 was applied by Affirmer to the Work. Should any part of the
License for any reason be judged legally invalid or ineffective under
applicable law, such partial invalidity or ineffectiveness shall not
invalidate the remainder of the License, and in such case Affirmer hereby
affirms that he or she will not (i) exercise any of his or her remaining
Copyright and Related Rights in the Work or (ii) assert any associated claims
and causes of action with respect to the Work, in either case contrary to
Affirmer's express Statement of Purpose.

4. Limitations and Disclaimers.

  a. No trademark or patent rights held by Affirmer are waived, abandoned,
  surrendered, licensed or otherwise affected by this document.

  b. Affirmer offers the Work as-is and makes no representations or warranties
  of any kind concerning the Work, express, implied, statutory or otherwise,
  including without limitation warranties of title, merchantability, fitness
  for a particular purpose, non infringement, or the absence of latent or
  other defects, accuracy, or the present or absence of errors, whether or not
  discoverable, all to the greatest extent permissible under applicable law.

  c. Affirmer disclaims responsibility for clearing rights of other persons
  that may apply to the Work or any use thereof, including without limitation
  any person's Copyright and Related Rights in the Work. Further, Affirmer
  disclaims responsibility for obtaining any necessary consents, permissions
  or other rights required for any use of the Work.

  d. Affirmer understands and acknowledges that Creative Commons is not a
  party to this document and has no duty or obligation with respect to this
  CC0 or use of the Work.

For more information, please see
<http://creativecommons.org/publicdomain/zero/1.0/>


text:
#!/bin/bash

# ./generate-interrupt-handler.sh > interrupt.S

for i in {0..255}
do
echo "
.global handle_interrupt_$i
.align 32
handle_interrupt_$i:
push rdi
mov rdi, $i
jmp handle_interrupt_"
done


text:
CFLAGS=-std=gnu99 -static -m64 -fno-PIC -nostdlib -nodefaultlibs -fno-stack-protector -U_FORTIFY_SOURCE -D_FORTIFY_SOURCE=0 -masm=intel -O1 -g -c

#.PHONY: musl

default: dep/musl/lib/libc.a dep/lua-5.2.3/src/liblua.a libsqlite3.a liblsqlite3.a
	./musl-custom-gcc $(CFLAGS) src/init.S -o bin/initS.o
	@# not sure why this has to be gcc instead of musl-gcc
	@# trap on lgdt if musl-gcc
	@# perhaps this is because musl is 64 bit but I pass -m32
	gcc $(CFLAGS) -m32 src/init.c -o bin/init.o
	@# init code is executed in protected mode, but is linked with 64-bit code
	objcopy -I elf32-i386 -O elf64-x86-64 bin/init.o bin/init.o
	@#~ luac -o luakernel.luac luakernel.lua
	@#~ xxd -i luakernel.luac luakernel.luac.h
	cd src && ../generate-lua-bundle.sh
	./musl-custom-gcc $(CFLAGS) -Idep/lua-5.2.3/src -Idep/sqlite3 src/luakernel.c -fno-stack-protector -o bin/luakernel.o
	./musl-custom-gcc -L. -T link.ld -fno-PIC -static -z max-page-size=0x1000 \
	  -o bin/luakernel.elf bin/initS.o bin/init.o bin/luakernel.o \
	  dep/lua-5.2.3/src/liblua.a libsqlite3.a liblsqlite3.a
	# make bootable ISO using GRUB
	rm -f bin/luakernel.iso bin/*.o
	@# comment out when debugging with gdb -- uncomment to get under 1.44 MB
# 	strip bin/luakernel.elf
	mkdir -p bin/boot/grub/
	cp grub.cfg bin/boot/grub
	@#pkgdatadir='~/scratch/grub/grub-core' ~/scratch/grub/grub-mkrescue -d /home/ers/scratch/grub/grub-core --locale-dir=/usr/lib/locale -o bin/luakernel.iso bin #-v

	grub-mkimage --format=i386-pc-eltorito --prefix="(cd)/boot/grub" --output=bin/grub.img \
	  --compression=none \
	  --config="bin/boot/grub/grub.cfg" loadenv biosdisk part_msdos part_gpt fat \
	  iso9660 loopback search boot minicmd cat cpuid chain \
	  halt help ls reboot echo test configfile normal sleep memdisk font \
	  gfxterm gettext true vbe vga video_bochs video_cirrus multiboot multiboot2

	genisoimage -graft-points -input-charset utf8 -A "LuaKernel" -quiet -R -b \
	  boot/grub/grub.img -no-emul-boot -boot-load-size 4 -boot-info-table \
	  -o bin/luakernel.iso \
	  luakernel.elf=bin/luakernel.elf \
	  boot/grub/grub.cfg=bin/boot/grub/grub.cfg \
	  boot/grub/grub.img=bin/grub.img

dep/musl/lib/libc.a:
	cd dep/musl && ./configure --disable-shared --enable-debug && make -j5

dep/lua-5.2.3/src/liblua.a:
	cd dep/lua-5.2.3 && make -j5 LDFLAGS="-m64 -fno-PIC" CC="gcc -m64 -O0 -g -fno-stack-protector -U_FORTIFY_SOURCE -D_FORTIFY_SOURCE=0 -fno-PIC -static" ansi

libsqlite3.a:
	gcc $(CFLAGS) -static -DSQLITE_ENABLE_MEMSYS5=1 -DSQLITE_THREADSAFE=0 dep/sqlite3/sqlite3.c -o libsqlite3.a

liblsqlite3.a:
	gcc -c $(CFLAGS) -Idep/lua-5.2.3/src/ -Idep/sqlite3 -static dep/lsqlite3.c -o liblsqlite3.a

run: default
	kvm -m 1024 -boot d -vga vmware -cdrom bin/luakernel.iso
	#qemu-system-x86_64 -kernel bin/luakernel.elf
	#~/scratch/bochs-2.6.6/bochs -f bochsrc.txt -q

dump:
	objdump -S -M intel -d bin/luakernel.elf | less

dumpdata:
	objdump -M intel -D bin/luakernel.elf | less

debug: default
	#kvm -s -m 1024 -boot d -vga vmware -cdrom bin/luakernel.iso
	#qemu-system-x86_64 -boot d -cdrom bin/luakernel.iso -s -S
	bochs-gdb-debugger -f bochsrc.txt -q
	#bochs-gui-debugger -f bochsrc.txt -q

gdb:
	#gdb -ex "set remote target-features-packet on" -ex "target remote localhost:1234" \
	#-ex "set architecture i386:x86-64:intel" bin/luakernel.elf

	#-ex "break trap" -ex "break resume_error" -ex "break luaG_runerror"

	gdb -ex "target remote localhost:1234" -ex "set remote target-features-packet on" \
	-ex "set architecture i386:x86-64:intel" \
	-ex "set disassembly-flavor intel" -ex "layout split" -ex "set confirm off" \
	-ex "break trap" \
	bin/luakernel.elf

	# p g->strt.size
	# p g->strt.hash[h % g->strt.size]
	# p *g->strt.hash@32
	# dump 32 quad words of string hash table
	# x /32xg 0x623710

	#-ex 'break internshrstr if str == 0x11b039' \

	#-ex 'break internshrstr if strcmp(str, "__gc") == 0' \

	# layout src
	# layout reg
	# layout asm
	# layout split

	# ctrl+x o to switch window focus for scrolling up in the command history

	#nemiver --remote=localhost:1234 --gdb-binary=/usr/bin/gdb bin/luakernel.elf
	#ddd --eval-command="target remote localhost:1234"

clean:
	rm -f bin/grub.img bin/luakernel.elf
	cd dep/musl && make clean
	cd dep/lua-5.2.3 && make clean


text:
# make dist-clean
# builtin graphical debugger
# ./configure --enable-debugger --enable-debugger-gui --enable-disasm --with-x \
# --enable-x86-64
# GDB
# ./configure --with-x --enable-x86-64 --enable-gdb-stub
# --with-nogui
megs: 1024
romimage: file=/usr/share/bochs/BIOS-bochs-latest
vgaromimage: file=/usr/share/bochs/VGABIOS-lgpl-latest
ata0-master: type=cdrom, path=bin/luakernel.iso, status=inserted
#display_library: sdl
#display_library: x, options="gui_debug"
display_library: x
#display_library: nogui
gdbstub: enabled=1, port=1234, text_base=0, data_base=0, bss_base=0
log: /tmp/out.bochs
mouse: enabled=1
cpu: ips=15000000, reset_on_triple_fault=0
boot: cdrom
magic_break: enabled=1
# http://bochs.sourceforge.net/doc/docbook/user/howto.html#LOGOPTS-BY-DEVICE
#debug: action=ignore, PIT=report
debug: action=ignore, cpu0=report
#debug: action=ignore, keyboard=report
#info: action=report, keyboard=ignore


text:
set timeout=0
set default=0
# set debug=all
menuentry "LuaKernel" {
	multiboot2 /luakernel.elf
	set gfxpayload=1280x768x32,1024x768x32,800x600x32,640x480x32
}


text:
#!/usr/bin/lua5.2

local raw = io.popen("find *.lua"):read("*all")
local filenames = {}
local modnames = {}
for filename in raw:gmatch("%a+%.lua") do
  table.insert(filenames, filename)
  table.insert(modnames, filename:sub(0, -5))
end
os.execute("rm -f lua-bundle.h")
for i, filename in ipairs(filenames) do
  os.execute(("xxd -i %s - >> lua-bundle.h"):format(filename))
end
local modname_str = ""
for i, modname in ipairs(modnames) do
  modname_str =
    modname_str .. ("\t{\"%s\", %s_lua, sizeof(%s_lua)},\n"):format(modname, modname, modname)
end
local array = [[
struct module
{
  char *name;
  unsigned char *buf;
  unsigned int len;
} lua_bundle[] =
{
%s
};
]]
local bundle = io.open("lua-bundle.h", "a+")
bundle:write(array:format(modname_str))


text:
#ifndef UTIL_H_
#define UTIL_H_

#include <stdint.h>

typedef int8_t s8;
typedef uint8_t u8;
typedef int16_t s16;
typedef uint16_t u16;
typedef int32_t s32;
typedef uint32_t u32;
typedef int64_t s64;
typedef uint64_t u64;

#define arraylen(array) (sizeof(array) / sizeof(array[0]))

void
outb(u32 addr, u32 value);

u32
inb(u32 addr);

#endif


text:
#include <fcntl.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <sys/syscall.h>
#include <sys/uio.h>
#include <time.h>
#include <unistd.h>
#include "lua.h"
#include "lualib.h"
#include "lauxlib.h"
#include "lua-bundle.h"
#include "multiboot2.h"
#include "sqlite3.h"
#include "util.h"
#include "vbe.h"

extern const uintptr_t heap_start;
static u8 *heap_end = 0;

static lua_State *L = NULL;

#if 1
// QEMU
void __attribute__((noinline))
trap()
{
  //~ asm volatile("");
  //~ asm volatile("int 3");
  asm volatile("xchg bx, bx");
}
#else
// bochs
#define trap() asm volatile("xchg bx, bx")
#endif

static void* l_alloc(void *ud, void *ptr, size_t osize, size_t nsize)
{
  (void)ud;
  (void)osize;
  if (nsize == 0)
  {
    sqlite3_free(ptr);
    return NULL;
  }
  else
  {
    return sqlite3_realloc(ptr, nsize);
  }
}

long handle_syscall(long n, long a1, long a2, long a3, long a4, long a5, long a6)
{
  switch (n)
  {
    case SYS_brk:
    {
      void *end = (void*)a1;
      // trap();
      return 0;
      break;
    }

    case SYS_open:
    {
      const char *pathname = (const char*)a1;
      int flags = a2;
      int mode = a3;

      lua_getglobal(L, "open");
      lua_pushstring(L, pathname);
      lua_pushnumber(L, flags);
      lua_pushnumber(L, mode);
      lua_call(L, 3, 1);
      return luaL_ref(L, LUA_REGISTRYINDEX);
      //int fd = luaL_ref(L, LUA_REGISTRYINDEX);
      //return lua_tonumber(L, 1);
      break;
    }

    case SYS_close:
    {
      return 0;
      break;
    }

    case SYS_lseek:
    {
      int fd = a1;
      off_t offset = a2;
      int whence = a3;

      lua_getglobal(L, "lseek");
      lua_rawgeti(L, LUA_REGISTRYINDEX, fd);
      lua_pushnumber(L, offset);
      lua_pushnumber(L, whence);
      lua_call(L, 3, 1);
      int retn = lua_tonumber(L, -1);
      lua_pop(L, 1);
      return retn;
      break;
    }

    case SYS_read:
    {
      int fd = a1;
      char *buf = (char*)a2;
      size_t count = a3;

      lua_getglobal(L, "read");
      if (fd != 0)
      {
        lua_rawgeti(L, LUA_REGISTRYINDEX, fd);
      }
      else
      {
        // stdin
        lua_pushnumber(L, fd);
      }
      lua_pushnumber(L, count);
      lua_call(L, 2, 1);
      size_t len;
      const char *resbuf = lua_tolstring(L, -1, &len);
      memcpy(buf, resbuf, len);
      lua_pop(L, 1);
      return len;
      //~ break;
    }

    case SYS_write:
    {
      int fd = a1;
      const char *buf = (const char*)a2;
      size_t count = a3;

      lua_getglobal(L, "write");
      if (fd != 1)
      {
        lua_rawgeti(L, LUA_REGISTRYINDEX, fd);
      }
      else
      {
        lua_pushnumber(L, fd);
      }
      lua_pushlstring(L, buf, count);
      lua_call(L, 2, 1);
      int len = lua_tonumber(L, -1);
      lua_pop(L, 1);
      return len;
    }

    case SYS_writev:
    {
      // p (char*)iov[0].iov_base
      int fd = a1;
      struct iovec *iov = (struct iovec*)a2;
      int count = (int)a3;
      ssize_t num_bytes_written = 0;
      for (int i = 0; i < count; ++i)
      {
        num_bytes_written += write(fd, iov[i].iov_base, iov[i].iov_len);
      }
      return num_bytes_written;
      break;
    }

    case SYS_ioctl:
    {
      return 0;
    }

    case SYS_getpid:
    {
      return 0;
    }

    case SYS_fsync:
    {
      return 0;
    }

    case SYS_chown:
    {
      return 0;
    }

    case SYS_clock_gettime:
    {
      clockid_t clk_id = (clockid_t)a1;
      struct timespec *tp = (struct timespec*)a2;
      *tp = (struct timespec){0};
      return 0;
    }

    case SYS_gettimeofday:
    {
      struct timeval *tv = (struct timeval*)a1;
      *tv = (struct timeval){0};
      return 0;
    }

    case SYS_fcntl:
    {
      return 0;
    }

    default:
    {
      trap();
      return -1;
    }
  }
}

long __syscall(long n, long a1, long a2, long a3, long a4, long a5, long a6)
{
  handle_syscall(n, a1, a2, a3, a4, a5, a6);
}

extern u8 *volatile multiboot_boot_information;
static struct VBEModeInfoBlock modeinfo;

static u8 *fbmem = NULL;
static u8 *display_buffer = NULL;
static u32 display_buffer_len = 0;

static int putpixel(lua_State *l)
{
  u32 x = lua_tonumber(l, 1);
  u32 y = lua_tonumber(l, 2);
  u32 r = lua_tonumber(l, 3);
  u32 g = lua_tonumber(l, 4);
  u32 b = lua_tonumber(l, 5);
  lua_pop(l, 5);
  const u32 bytes_per_pixel = (modeinfo.BitsPerPixel / 8);
  // http://forum.osdev.org/viewtopic.php?p=77998&sid=d4699cf03655c572906144641a98e4aa#p77998
  u8 *ptr =
    &display_buffer[(y * modeinfo.BytesPerScanLine) + (x * bytes_per_pixel)];
  const u8 *display_buffer_end =
    &display_buffer[(modeinfo.YResolution * modeinfo.BytesPerScanLine)];
  if (ptr < display_buffer_end)
  {
    ptr[0] = b;
    ptr[1] = g;
    ptr[2] = r;
    ptr[3] = 0;
  }
  else
  {
    trap();
  }

  return 0;
}

static int clear_screen(lua_State *l)
{
  memset(display_buffer, 0, display_buffer_len);
  return 0;
}

static int swap_buffers(lua_State *l)
{
  memcpy(fbmem, display_buffer, display_buffer_len);
  clear_screen(l);
  return 0;
}

static void get_multiboot_info(void)
{
  if (multiboot_boot_information == NULL)
  {
    trap();
  }
  u32 total_size = *(u32*)multiboot_boot_information;
  struct multiboot_tag *tag = (struct multiboot_tag*)&multiboot_boot_information[8];
  while (tag->type != MULTIBOOT_TAG_TYPE_END)
  {
    switch (tag->type)
    {
      case MULTIBOOT_TAG_TYPE_VBE:
      {
        struct multiboot_tag_vbe *vbetag = (struct multiboot_tag_vbe*)tag;
        modeinfo = *(struct VBEModeInfoBlock*)&vbetag->vbe_mode_info;
        break;
      }

      case MULTIBOOT_TAG_TYPE_FRAMEBUFFER:
      {
        struct multiboot_tag_framebuffer *fb = (struct multiboot_tag_framebuffer*)tag;
        if (modeinfo.XResolution > 0)
        {
          fbmem = (u8*)fb->common.framebuffer_addr;
        }
        break;
      }
#if 0
      case MULTIBOOT_TAG_TYPE_BASIC_MEMINFO:
      {
        struct multiboot_tag_basic_meminfo *meminfo = (struct multiboot_tag_basic_meminfo*)tag;
        break;
      }
#endif

      case MULTIBOOT_TAG_TYPE_MMAP:
      {
        struct multiboot_tag_mmap *mmap = (struct multiboot_tag_mmap*)tag;
        for
        (
          struct multiboot_mmap_entry *entry = mmap->entries;
          (u8*)entry < (u8*)mmap + tag->size;
          entry = (struct multiboot_mmap_entry*)((u8*)entry + mmap->entry_size)
        )
        {
          if (entry->type == MULTIBOOT_MEMORY_AVAILABLE)
          {
            // available memory:
            // 0x0 to 0x9f000
            // 0x100000 to 0x3FFF0000
            // 0xdcaa00
            if (entry->addr == 0x100000)
            {
              heap_end = (u8*)entry->addr + entry->len;
              //~ trap();
            }
          }
        }
        break;
      }
    }

    // tags are padded to ensure 8 byte alignment
    if (tag->size % 8 == 0)
    {
      tag = (struct multiboot_tag*)((u8*)tag + tag->size);
    }
    else
    {
      tag = (struct multiboot_tag*)((u8*)tag + tag->size + (8 - (tag->size % 8)));
    }
  }
}

static u64 timer_ticks = 0;

static u32 keyboard_scancode_queue[8] = {0};
static u32 keyboard_scancode_queue_len = 0;

void handle_interrupt(u32 n)
{
  switch (n)
  {
    // general protection fault
    case 13:
    {
      trap();
      break;
    }

    // page fault
    case 14:
    {
      trap();
      break;
    }

    // timer
    case 32:
    {
      //~ trap(); while (1);
      // 100 Hz
      timer_ticks += 1;
      outb(0x20, 0x20);
      break;
    }

    // keyboard
    case 33:
    {
      u32 scancode = inb(0x60);
      if (keyboard_scancode_queue_len < arraylen(keyboard_scancode_queue))
      {
        keyboard_scancode_queue[keyboard_scancode_queue_len] = scancode;
        keyboard_scancode_queue_len += 1;
      }
      outb(0x20, 0x20);
      break;
    }

    // mouse
    case 44:
    {
      u32 n = inb(0x60);
      outb(0xa0, 0x20);
      outb(0x20, 0x20);
      break;
    }

    default:
    {
      trap(); while (1);
      break;
    }
  }
}

static int lua_outb(lua_State *l)
{
  const u32 addr = lua_tonumber(l, 1);
  const u8 value = lua_tonumber(l, 2);
  outb(addr, value);
  lua_pop(l, 2);
  return 0;
}

static int lua_inb(lua_State *l)
{
  const u32 addr = lua_tonumber(l, 1);
  lua_pop(l, 1);
  const u8 value = inb(addr);
  lua_pushnumber(l, value);
  return 1;
}

// http://lua-users.org/lists/lua-l/2003-12/msg00301.html
// http://lua-users.org/lists/lua-l/2002-12/msg00171.html
// http://lua-users.org/lists/lua-l/2011-06/msg00426.html
// http://lua-users.org/lists/lua-l/2010-03/msg00679.html

static void lua_hook(lua_State *l, lua_Debug *ar)
{
  lua_yield(l, 0);
}

static int lua_setmaskhook(lua_State *l)
{
  lua_State *t = lua_tothread(l, 1);
  int maskcount = lua_tointeger(l, 2);
  lua_pop(l, 2);
  if (t)
  {
    lua_sethook(t, lua_hook, LUA_MASKCOUNT, maskcount);
  }
  return 0;
}

static int lua_get_timer_ticks(lua_State *l)
{
  lua_pushinteger(l, timer_ticks);
  return 1;
}

static int lua_get_keyboard_interrupt(lua_State *l)
{
  // disable interrupts
  asm volatile ("cli");

  // process interrupt data
  lua_createtable(l, keyboard_scancode_queue_len, 0);
  for (int i = 0; i < keyboard_scancode_queue_len; ++i)
  {
    lua_pushinteger(l, keyboard_scancode_queue[i]);
    lua_rawseti(l, -2, i + 1);
  }
  keyboard_scancode_queue_len = 0;

  // enable interrupts
  asm volatile ("sti");
  return 1;
}

static int lua_hlt(lua_State *l)
{
  asm volatile("hlt");
  return 0;
}

const char *errstr = NULL;

static int lua_loader(lua_State *l)
{
  size_t len;
  const char *modname = lua_tolstring(l, -1, &len);
  struct module *mod = NULL;
  for (int i = 0; i < arraylen(lua_bundle); ++i)
  {
    if (memcmp(modname, lua_bundle[i].name, len) == 0)
    {
      mod = &lua_bundle[i];
    }
  }
  if (!mod)
  {
    lua_pushnil(l);
    return 1;
  }
  if (luaL_loadbuffer(l, mod->buf, mod->len, mod->name) != LUA_OK)
  {
    errstr = lua_tostring(l, 1);
    //~ puts("luaL_loadstring: error");
    trap();
  }
  int err = lua_pcall(l, 0, LUA_MULTRET, 0);
  if (err != LUA_OK)
  {
    errstr = lua_tostring(l, 1);
    //~ puts("lua_pcall: error");
    trap();
  }
  if (!lua_istable(l, -1))
  {
    puts("not a table");
  }
  return 1;
}

void main(void)
{
  get_multiboot_info();

  // Use SQLite3 as the only memory allocator because musl's malloc requires mmap.
  sqlite3_config(SQLITE_CONFIG_HEAP, &heap_start, heap_end - (u8*)&heap_start, 64);

  L = lua_newstate(l_alloc, NULL);
  if (!L)
  {
    puts("lua_newstate: error");
    return;
  }
  luaL_openlibs(L);

  display_buffer_len = (modeinfo.YResolution * modeinfo.BytesPerScanLine);
  display_buffer = lua_newuserdata(L, display_buffer_len);
  clear_screen(L);

  lua_pushnumber(L, modeinfo.XResolution);
  lua_setglobal(L, "DISPLAY_WIDTH");
  lua_pushnumber(L, modeinfo.YResolution);
  lua_setglobal(L, "DISPLAY_HEIGHT");
  lua_register(L, "clear_screen", clear_screen);
  lua_register(L, "putpixel", putpixel);
  lua_register(L, "swap_buffers", swap_buffers);
  lua_register(L, "outb", lua_outb);
  lua_register(L, "inb", lua_inb);
  lua_register(L, "setmaskhook", lua_setmaskhook);
  lua_register(L, "loader", lua_loader);
  lua_register(L, "get_timer_ticks", lua_get_timer_ticks);
  lua_register(L, "get_keyboard_interrupt", lua_get_keyboard_interrupt);
  lua_register(L, "hlt", lua_hlt);

  int luaopen_lsqlite3(lua_State *L);
  luaL_requiref(L, "lsqlite3", luaopen_lsqlite3, 0);
  lua_pop(L, 1);

  if (luaL_loadbuffer(L, luakernel_lua, luakernel_lua_len, "luakernel") != LUA_OK)
  {
    //~ puts("luaL_loadstring: error");
    errstr = lua_tostring(L, 1);
    trap();
    return;
  }
  int err = lua_pcall(L, 0, LUA_MULTRET, 0);
  if (err != LUA_OK)
  {
    //~ puts("lua_pcall: error");
    trap();
    return;
  }
  trap();
}


text:
// https://github.com/mokafive/grub/blob/upstream/grub-core/video/i386/pc/vbe.c

struct VBEModeInfoBlock
{
  // Mandatory information for all VBE revisions
  u16 ModeAttributes;
  u8 WinAAttributes;
  u8 WinBAttributes;
  u16 WinGranularity;
  u16 WinSize;
  u16 WinASegment;
  u16 WinBSegment;
  u32 WinFuncPtr;
  u16 BytesPerScanLine;
  // Mandatory information for VBE 1.2 and above
  u16 XResolution;
  u16 YResolution;
  u8 XCharSize;
  u8 YCharSize;
  u8 NumberOfPlanes;
  u8 BitsPerPixel;
  u8 NumberOfBanks;
  u8 MemoryModel;
  u8 BankSize;
  u8 NumberOfImagePages;
  u8 Reserved_page;
  // Direct Color fields (required for direct/6 and YUV/7 memory models)
  u8 RedMaskSize;
  u8 RedFieldPosition;
  u8 GreenMaskSize;
  u8 GreenFieldPosition;
  u8 BlueMaskSize;
  u8 BlueFieldPosition;
  u8 RsvdMaskSize;
  u8 RsvdFieldPosition;
  u8 DirectColorModeInfo;
  // Mandatory information for VBE 2.0 and above
  u32 PhysBasePtr;
  u32 OffScreenMemOffset;
  u16 OffScreenMemSize;
  // Mandatory information for VBE 3.0 and above
  u16 LinBytesPerScanLine;
  u8 BnkNumberOfPages;
  u8 LinNumberOfPages;
  u8 LinRedMaskSize;
  u8 LinRedFieldPosition;
  u8 LinGreenMaskSize;
  u8 LinGreenFieldPosition;
  u8 LinBlueMaskSize;
  u8 LinBlueFieldPosition;
  u8 LinRsvdMaskSize;
  u8 LinRsvdFieldPosition;
  u32 MaxPixelClock;
  u8 Reserved[189];
} __attribute__((packed));