using GaussianProcesses
using Random

Random.seed!(20140430)
# Training data
n=10;                          #number of training points
x = 2π * rand(n);              #predictors
y = sin.(x) + 0.05*randn(n);   #regressors

#Select mean and covariance function
mZero = MeanZero()                   #Zero mean function
kern = SE(0.0,0.0)                   #Squared exponential kernel (note that hyperparameters are on the log scale)

logObsNoise = -1.0                        # log standard deviation of observation noise (this is optional)
gp = GP(x,y,mZero,kern,logObsNoise)       #Fit the GP

μ, σ² = predict_y(gp,range(0,stop=2π,length=100));

using Plots  #Load Plots.jl package

plot(gp; xlabel="x", ylabel="y", title="Gaussian process", legend=false, fmt=:png)      # Plot the GP

using Optim
optimize!(gp; method=ConjugateGradient())   # Optimise the hyperparameters

plot(gp; legend=false, fmt=:png)   #Plot the GP after the hyperparameters have been optimised 

optimize!(gp; kern = false)   # Don't optimize kernel hyperparameters
optimize!(gp; kernbounds = [[-1, -1], [1, 1]]) # Optimize the kernel parameters in a box with lower bounds [-1, -1] and upper bounds [1, 1]

using Distributions

set_priors!(kern, [Normal(), Normal()]) # Uniform(0,1) distribution assumed by default if priors not specified
@time chain = mcmc(gp, nIter=5000);

set_priors!(gp.logNoise, [Normal()])
@time ess_chain = ess(gp, nIter=5000);

p1 = plot(chain', label=["Noise" "SE log length" "SE log scale"]; fmt=:png, size = (1200, 600), title="Hamiltonian Monte-Carlo", ylims=(-4, 2.5))
p2 = plot(ess_chain', label=["Noise" "SE log length" "SE log scale"]; fmt=:png, size = (1200, 600), title="Elliptical Slice Sampler", ylims=(-4, 2.5))
plot(p1, p2, xlabel="Iteration", ylabel="Parameter Value")

#Training data
d, n = 2, 50;         #Dimension and number of observations
x = 2π * rand(d, n);                               #Predictors
y = vec(sin.(x[1,:]).*sin.(x[2,:])) + 0.05*rand(n);  #Responses

mZero = MeanZero()                             # Zero mean function
kern = Matern(5/2,[0.0,0.0],0.0) + SE(0.0,0.0)    # Sum kernel with Matern 5/2 ARD kernel 
                                               # with parameters [log(ℓ₁), log(ℓ₂)] = [0,0] and log(σ) = 0
                                               # and Squared Exponential Iso kernel with
                                               # parameters log(ℓ) = 0 and log(σ) = 0

gp = GP(x,y,mZero,kern,-2.0)          # Fit the GP

optimize!(gp)                         # Optimize the hyperparameters

plot(contour(gp) ,heatmap(gp); fmt=:png)