library(MASS)

l=1
sigma_f=1

k1<-function(x1,x2){return(sigma_f*exp(-0.5*sum((x1-x2)^2)/l^2))}
k<-function(x1,x2){return((x1*x2+2)^2)}
k3<-function(x1,x2){return(exp(-abs(x1-x2)))}

xPoints<-c(-1,0,1)
y<-c(1,0,1)


Kernel = matrix( , nrow=length(xPoints), ncol=length(xPoints) )

for (i in 1:length(xPoints) ) {
    for (j in 1:length(xPoints) ) {
      Kernel[i,j] = k(xPoints[i],xPoints[j]) 
#       Kernel[j,i] = Kernel[i,j]
    }
}
KInverse = solve(Kernel)

predicted=c()
predictedPlusVariance=c()
predictedMinusVariance=c()

xSeq = seq(-5,5,by=0.1)

groundTruth=xSeq^2

for (i in 1:length(xSeq)) {
  kStar=c()
  for(j in 1:length(xPoints) ) {
    kStar = c(kStar,k(xSeq[i],xPoints[j]))
  }
#   kStar = c(k(xSeq[i],xPoints[1]),k(xSeq[i],xPoints[2]))
  kStarStar = k(xSeq[i],xSeq[i])
  mean = t(kStar)%*%KInverse%*%y
  predicted = c(predicted,mean)
  var = kStarStar - t(kStar)%*%KInverse%*%kStar
  predictedPlusVariance = c(predictedPlusVariance,mean+var)
  predictedMinusVariance = c(predictedMinusVariance,mean-var)
}

plot(xSeq,ylim=c(-0.5,2),predicted,type="l")
lines(xSeq,predictedMinusVariance,type="l",col="blue")
lines(xSeq,predictedPlusVariance,type="l",col="blue")
lines(xSeq,groundTruth,type="l",col="red")


# 
# updateKernelMatrix<-function(x)
# {
# 	K<-matrix(ncol=nrow(x))
# 	for(i in 1:nrow(x))
# 	{
# 		K<-rbind(K,apply(x,1,k,x2=x[i]))
# 	}
# 	return(K[-1,])
# }
# 
# gpMean<-function(alpha,xStar,x)
# {
# 	kStar<-apply(x,1,k,x2=xStar)
# 	mean<-kStar%*%alpha
# 	return(mean)
# }
# 
# gpSd<-function(K,xStar,x)
# {
# 	kStarStar<-k(xStar,xStar)
# 	kStar<-apply(x,1,k,x2=xStar)
# 	sd<-kStarStar-kStar%*%solve(K)%*%kStar
# 	return(sd)
# }
# 
# 
# K<-updateKernelMatrix(x)
# alpha<-solve(K)%*%y
# 
# plot(x,y,xlim=c(-10,10),ylim=c(-1,1.5))
# 
# a<-matrix(seq(-10,10,0.1),ncol=1)
# lines(a,apply(a,1,gpMean,x=x,alpha=alpha),col="red")
# 
# lines(a,apply(a,1,gpMean,x=x,alpha=alpha)+apply(a,1,gpSd,x=x,K=K),col="gray",lty=2)
# lines(a,apply(a,1,gpMean,x=x,alpha=alpha)-apply(a,1,gpSd,x=x,K=K),col="gray",lty=2)
# 
# x<-rbind(x,0)
# y<-rbind(y,0)
# K<-updateKernelMatrix(x)
# alpha<-solve(K)%*%y
# dev.new()
# plot(x,y,xlim=c(-10,10),ylim=c(-1,1.5))
# 
# a<-matrix(seq(-10,10,0.01),ncol=1)
# lines(a,apply(a,1,gpMean,x=x,alpha=alpha),col="red")
# 
# lines(a,apply(a,1,gpMean,x=x,alpha=alpha)+apply(a,1,gpSd,x=x,K=K),col="gray",lty=2)
# lines(a,apply(a,1,gpMean,x=x,alpha=alpha)-apply(a,1,gpSd,x=x,K=K),col="gray",lty=2)



# Exercise 3:

# 
# xPoints=seq(0,5,by=0.1)
# mean = rep(0,length(xPoints))

# Kernel = matrix( , nrow=length(xPoints), ncol=length(xPoints) )
# 
# for (i in 1:length(xPoints) ) {
#     for (j in i:length(xPoints) ) {
#       Kernel[i,j] = k1(xPoints[i],xPoints[j])
#       Kernel[j,i] = Kernel[i,j]
#     }
# }
# 
# f = mvrnorm(1,mean,Kernel)
# plot(xPoints,f,type="l",col=1)
# 
# for ( i in 1:5) {
#   f = mvrnorm(1,mean,Kernel)
#   lines(xPoints,f,col=i+1)
# }