library(pracma)
library(mvtnorm)
library(rcage)



covar = function(x, y, tau=1){
  return( exp(-1/tau *abs(x-y)))
}
n= 200
rho = 1
x = seq(-1, 1, length.out = n)



compute_KLE = function(nbasis, rho){
  
  P = list()
  P[[1]]=1;
  P[[2]]=c(1, 0)
  for (i in 2:nbasis){#recurrence relation{
    P[[i+1]] =1/i* c( c((2*i-1)*P[[i]], 0) - (i-1)*c(0, 0, P[[i-1]])    )
  }
  
  for (i in 1:nbasis){#normalization
    P[[i]] = P[[i]] / sqrt(2/(2*i-1));
  } 
  
  
  
  
  Theta=zeros(n, nbasis)
  for (i in 1:nbasis){
    Theta[,i] = polyval(P[[i]],x)
  }
  
  A=zeros(nbasis);
  B=zeros(nbasis);
  
  
  
  for (i in 1:nbasis){
    for (j in 1:nbasis){
      B[i,j]=trapz(x, Theta[,i] * Theta[,j]);
    }
  }
  
  
  for (i in 1:nbasis){
    print(i)
    for (j in 1:nbasis){
      temp = c()
      for (s in 1:n) {
        temp = c(temp, trapz(x, covar(x[s], x, rho) * Theta[,j]))
      }
      A[i,j]=trapz(x, temp * Theta[, i])
    }
  }
  
  
  G = solve(B) %*% A
  eig = eigen(G)
  lam = eig$values
  F = eig$vectors
  Psi = Theta %*% F
  
  for (i in 1:5) {
    for (j in 1:5) {
      cat(i, j, trapz(x, Psi[,i] * Psi[,j]), "\n")
    }
  }
  
  Lam = diag(lam)
  H = Psi %*% Lam %*% t(Psi)
  XY = meshgrid(x,x)
  C = covar(XY$X, XY$Y, rho)
  return(list(Psi=Psi, Lam=Lam, H=H, C=C))
}


KLE1_10 = compute_KLE(10, rho)
KLE1_5 = compute_KLE(5, rho)
KLE1_15 = compute_KLE(15, rho)


XY = meshgrid(x,x)





png("legendre1d_ev.png")
plot(diag(KLE1_5$Lam), type="l", col="red", xlim=c(1, 10), ylim=c(0, 1.16), ylab = "", xlab = "", main = "Eigenvalues")
par(new=T)
plot(diag(KLE1_10$Lam), type="l", col="green", xlim=c(1, 10), ylim=c(0, 1.15), ylab = "", xlab = "", main = "Eigenvalues")
legend(x = "topright",          # Position
       legend = c("5 basis", "10 basis"),  # Legend texts
       col = c("red", "green"),           # Line colors
       lwd = 1)  
dev.off()


png("legendre1d_ef.png")
par(mfrow=c(1,2))
#plot(x, KLE1_5$Psi[,1], type = "l", main = "5 Basis", xlab = "", ylab = "First Eigenfunction")
#plot(x, KLE1_5$Psi[,2], type = "l", main = "5 Basis", xlab = "", ylab = "Second Eigenfunction")
plot(x, KLE1_10$Psi[,1], type = "l",  xlab = "x", ylab=expression(paste(psi[1], "(x)")), main = "First Eigenfunction")
plot(x, KLE1_10$Psi[,2], type = "l",  xlab = "x", ylab=expression(paste(psi[2], "(x)")),  main = "Second Eigenfunction")
dev.off()


