### Rewrite using makeLLS and useLLS conventions, and standardize

makeBayesReg <- function(X, T, alpha, beta) {
  X1 <- cbind(1,X)
  M <- ncol(X1)
  S.w <- solve(beta * t(X1) %*% X1 + alpha * diag(1,M))
  mu.w <- beta * S.w %*% t(X1) %*% T
  list(mu.w = mu.w, S.w = S.w, alpha = alpha, beta = beta)
}

useBayesReg <- function(model, X) {
  X1 <- cbind(1,X)
  mu.t <- X1 %*% model$mu.w
  S.t <- 1/model$beta + rowSums(X1 %*% model$S.w * X1)
browser()
  append(list(mean = mu.t, variance = S.t), model)
}

plot.BayesRegression1D <- function(model,Xtrain,Ttrain,Xtest,predictions,nRBFs,rbfWidth) {
  upper <- predictions$mean + sqrt(predictions$variance)
  lower <- predictions$mean - sqrt(predictions$variance)
  ## See demo(graphics) to see how to draw shaded region
  plot(c(Xtest,rev(Xtest)),c(upper,rev(lower)),type="n",xlab="x",ylab="t",
       ylim=c(-2,2))
  polygon(c(Xtest,rev(Xtest)),c(upper,rev(lower)),col="pink",border=NA)
  points(Xtrain,Ttrain,col="blue")
  lines(Xtest,predictions$mean,col="red")
  lines(Xtest,sin(2*pi*Xtest),col="green")
  title(paste("nRBFs =",nRBFs," rbfWidth =",rbfWidth,
              " nTrain =",length(Xtrain),"\n",
              " beta =",model$beta,
              " alpha = ",model$alpha))
}

makeStandardizeF <- function(X) {
  if (missing(X)) {
    cat("Usage:
         standardize <- makeStandardizeF(X)  ## X is nSamples x nDimensions
         Xs <- standardize(X)
         X2s <- standardize(X2)\n")
    return(invisible())
  }
  ## X is nSamples x nDimensions
  mu <- colMeans(X)
  sigma <- sd(X) ##sd should be named colSds

  function(newX) {
    nr <- nrow(newX)
    nc <- ncol(newX)
    (newX - matrix(mu,nr,nc,byrow=TRUE)) / matrix(sigma,nr,nc,byrow=TRUE)
  }
}

######################################################################
### Make the dta

rbfize1D <- function(X,centers,widths) {
  ## X is column matrix
  nSamples <- nrow(X)
  nbf <- length(centers)
  exp(-(matrix(X,nSamples,nbf) - matrix(centers,nSamples,nbf,byrow=TRUE))^2 /
      matrix(widths,nSamples,nbf,byrow=TRUE))
}

### First few picked for nTrain=1,2, or 4
Xall <- matrix(c(0.4,0.6,0,1,runif(50)))
nRBFs <- 9
rbfWidth <- 0.05
xrange <- range(Xall)
rbfCenters <- seq(xrange[1],xrange[2],len=nRBFs)

### Convert X to RBF values
XallRBF <- rbfize1D(Xall, centers=rbfCenters, widths=rbfWidth)
beta <- 25
Tall <- matrix(sin(2*pi*Xall)+rnorm(nrow(Xall),0,sqrt(1/beta)))  ## assumes beta is known!
nTest <- 100
Xtest <- matrix(seq(0,1,len=nTest))
Ttest <- matrix(sin(2*pi*Xtest)+rnorm(nrow(Xtest),0,sqrt(1/beta)))
XtestRBF <- rbfize1D(Xtest, centers=rbfCenters, widths=rbfWidth)

### Precision Parameters
alpha <- 0.1

### Set up graphics
def.par <- par(mar=c(2,2,3,1))
layout(matrix(c(1,2,3,4),2,2,byrow=TRUE))  ## read ?layout

### Use four different values of nTrain
for (nTrain in c(1,2,4,25)) {

  Xtrain <- Xall[1:nTrain,,drop=FALSE]
  XtrainRBF <- XallRBF[1:nTrain,,drop=FALSE]
  Ttrain <- Tall[1:nTrain,,drop=FALSE]

  model <- makeBayesReg(XtrainRBF,Ttrain,alpha,beta)
  predictions <- useBayesReg(model,XtestRBF)

  plot.BayesRegression1D(model,Xtrain,Ttrain,Xtest,predictions,nRBFs,rbfWidth)
}

par(def.par)

######################################################################
### now the auto mpg data

mpg <-read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data")
mpgnames <- c("mpg","cylinders","displacement","horsepower",
              "weight","acceleration","year","origin","name")
colnames(mpg) <- mpgnames

### Remove all samples that have at least one "?"
keepRows <- apply(mpg != "?", 1, all)
mpg <- mpg[keepRows,]

### Randomly pick 80% of samples for training partition
randorder <- sample(nrow(mpg))
nTrain <- round(nrow(mpg)*0.8)
trainRows <- randorder[1:nTrain]

### Assemble X and T matrices
Xtrain <- mpg[trainRows, 2:8]
Ttrain <- mpg[trainRows, 1, drop=FALSE]

Xtrain <- apply(Xtrain,2,as.numeric)
Ttrain <- apply(Ttrain,2,as.numeric)
standardize <- makeStandardizeF(Xtrain)
Xtrain <- standardize(Xtrain)

Xtest <- mpg[-trainRows, 2:8]
Ttest <- mpg[-trainRows, 1, drop=FALSE]
Xtest <- apply(Xtest,2,as.numeric)
Ttest <- apply(Ttest,2,as.numeric)
Xtest <- standardize(Xtest)

model <- makeBayesReg(Xtrain,Ttrain,alpha,beta)
predictions <- useBayesReg(model,Xtest)

for (i in 1:10) {
  cat(predictions$mean[i],"+-",sqrt(predictions$variance[i]),"\n") 
}