using Plots
pyplot()   # Optionally select a plotting backend

using GaussianProcesses

# Generate random data for Gaussian process
x = 2π * rand(10)
y = sin.(x) + 0.5*rand(10)

# Set-up mean and kernel
se = SE(0.0, 0.0)
m = MeanZero()

# Construct and plot GP
gp = GP(x,y,m,se)
plot(gp;  xlabel="gp.x", ylabel="gp.y", title="Gaussian process", legend=false, fmt=:png)

x = 0:0.1:2π
plot(gp; obsv=false)
optimize!(gp)
plot(gp; obsv=false, label="GP posterior mean", fmt=:png)
samples = rand(gp, x, 5)
plot!(x, samples)

# Simulate data for 2D Gaussian process
X = 2π*rand(2, 10)
y = sin.(X[1,:]) .* cos.(X[2,:]) + 0.5*rand(10)
gp2 = GP(X,y,m,se)
# Plot mean and variance
p1 = plot(gp2; title="Mean of GP")
p2 = plot(gp2; var=true, title="Variance of GP", fill=true)
plot(p1, p2; fmt=:png)

gr() # use GR backend to allow wireframe plot
p1 = contour(gp2)
p2 = surface(gp2)
p3 = heatmap(gp2)
p4 = wireframe(gp2)
plot(p1, p2, p3, p4; fmt=:png)