X <- matrix(seq(0,1,len=10))
#f <- function(X) matrix(sin(X) * sin(0.8*X/1.4 + 5) + rnorm(length(X),0,0.1))
f <- function(X) matrix(sin(3*X) + rnorm(length(X),0,0.2))
T <- f(X)
Xtest <- X +0.05
Ttest <- f(Xtest)

x11()
plot(X,T,type="b",main="Data",ylim=range(c(T,Ttest)),bty="n")
lines(Xtest,Ttest,type="b",col=2)
legend("topright",c("Train","Test"),lty=1,col=1:2)
dev.copy2pdf(file="noisysine.pdf")

polynomialize <- function(X,degrees=1:3) {
  Xnew <- c()
  for (d in degrees) {
    Xnew <- cbind(Xnew, X^d)
  }
}

rbfize <- 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))
}


x11()
xs <- seq(0,1,length=50)
phis <- polynomialize(xs,1:3)
matplot(xs,phis,type="l",xlab="x",ylab=expression(phi[d](x)),col=c(1,2,4),lwd=2)
dev.copy2pdf(file="polynomials.pdf")

phis <- rbfize(matrix(xs),seq(0,1,by=0.2),0.1)
matplot(xs,phis,type="l",xlab="x",ylab=expression(phi[d](x)),lwd=2)
dev.copy2pdf(file="rbfs.pdf")


x11()
p <- par(mfrow=c(3,3))
results <- c()
for (maxDegree in 1:9) {

  Xpoly <- polynomialize(X,1:maxDegree)
  Xtestpoly <- polynomialize(Xtest,1:maxDegree)

  X1 <- cbind(1,Xpoly)
  w <- solve(t(X1) %*% X1, t(X1) %*% T)

  X1test <- cbind(1,Xtestpoly)
  preds <- X1test %*% w

  plot(X,T,type="b",col=1,main=paste("Degree",maxDegree),bty="n")
  lines(Xtest,Ttest,type="b",col=2)
  lines(Xtest,preds,type="b",col=4)

  rmse <- sqrt(mean((preds-Ttest)^2))
  results <- rbind(results,
                   c(maxDegree,rmse))

}
dev.copy2pdf(file="multipleDegrees.pdf")

x11()
plot(results[,1],results[,2],type="b",xlab="Polynomial Degree",ylab="Test RMSE",,bty="n")
dev.copy2pdf(file="multipleDegreesResults.pdf")



x11()

p <- par(mfrow=c(3,3))
results <- c()
for (numRBFs in 1:8) {

  centers <- seq(0,1,len=numRBFs)
  widths <- rep(0.1,length(centers))
  Xpoly <- rbfize(X,centers,widths)
  Xtestpoly <- rbfize(Xtest,centers,widths)

  X1 <- cbind(1,Xpoly)
  w <- solve(t(X1) %*% X1, t(X1) %*% T)

  X1test <- cbind(1,Xtestpoly)
  preds <- X1test %*% w

  plot(X,T,type="b",col=1,main=paste("Num RBFs",numRBFs),bty="n")
  lines(Xtest,Ttest,type="b",col=2)
  lines(Xtest,preds,type="b",col=4)

  rmse <- sqrt(mean((preds-Ttest)^2))
  results <- rbind(results,
                   c(numRBFs,rmse))
}
dev.copy2pdf(file="multipleRBFs.pdf")

x11()
plot(results[,1],results[,2],type="b",xlab="Number RBFs",ylab="Test RMSE",bty="n")
dev.copy2pdf(file="multipleRBFsResults.pdf")