remove stheno; add AbstractGP surrogate

Project.toml

@@ -4,13 +4,17 @@ authors = ["ludoro <[email protected]>"]
 version = "4.1.0"
 
 [deps]
+AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
+Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LatinHypercubeSampling = "a5e1c1ea-c99a-51d3-a14d-a9a37257b02d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+PolyChaos = "8d666b04-775d-5f6e-b778-5ac7c70f65a3"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 XGBoost = "009559a3-9522-5dbb-924b-0b6ed2b22bb9"
@@ -26,7 +30,6 @@ LatinHypercubeSampling = "1.2"
 PolyChaos = "0.2"
 Requires = "0.5, 1.0"
 Sobol = "1.3"
-Stheno = "0.6"
 XGBoost = "0.4, 1.1"
 Zygote = "0.4, 0.5, 0.6"
 julia = "1.6"

docs/make.jl

@@ -7,6 +7,7 @@ makedocs(
         "Basics" => "tutorials.md",
         "Radials" => "radials.md",
         "Kriging" => "kriging.md",
+        "Gaussian Process" =>"abstractgps.md",
         "Lobachevsky" => "lobachevsky.md",
         "Linear" => "LinearSurrogate.md",
         "InverseDistance" => "InverseDistance.md",
@@ -17,7 +18,7 @@ makedocs(
         "Polynomial Chaos" => "polychaos.md",
         "Variable Fidelity" => "variablefidelity.md",
         "Gradient Enhanced Kriging" => "gek.md"
-        ]
+       ]
     "User guide" => [
         "Samples" => "samples.md",
         "Surrogates" => "surrogate.md",

docs/src/abstractgps.md

@@ -0,0 +1,100 @@
+# Gaussian Process Surrogate Tutorial
+
+Gaussian Process regression in Surrogates.jl is implemented as a simple wrapper around the [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) package. AbstractGPs comes with a variety of covariance functions (kernels). See [KernelFunctions.jl](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/) for examples.
+
+
+## 1D Example 
+In the example below, the 'gp_surrogate' assignment code can be commented / uncommented to see how the different kernels influence the predictions. 
+
+```@example gp_tutorial1d
+using Surrogates
+using Plots
+default()
+using AbstractGPs #required to access different types of kernels
+
+f(x) = (6 * x - 2)^2 * sin(12 * x - 4)
+n_samples = 4
+lower_bound = 0.0
+upper_bound = 1.0
+xs = lower_bound:0.001:upper_bound
+x = sample(n_samples, lower_bound, upper_bound, SobolSample())
+y = f.(x)
+#gp_surrogate = AbstractGPSurrogate(x,y, gp=GP(SqExponentialKernel()), Σy=0.05) #example of Squared Exponential Kernel
+#gp_surrogate = AbstractGPSurrogate(x,y, gp=GP(MaternKernel()), Σy=0.05) #example of MaternKernel
+gp_surrogate = AbstractGPSurrogate(x,y, gp=GP(PolynomialKernel(; c=2.0, degree=5)), Σy=0.25)
+plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound), ylims=(-7, 17), legend=:top)
+plot!(xs, f.(xs), label="True function", legend=:top)
+plot!(0:0.001:1, gp_surrogate.gp_posterior; label="Posterior", ribbon_scale=2)
+```
+
+## Optimization Example
+This example shows the use of AbstractGP Surrogates to find the minima of a function:
+
+```@example abstractgps_tutorial_optimization
+using Surrogates
+using Plots
+f(x) = (x-2)^2
+n_samples = 4
+lower_bound = 0.0
+upper_bound = 4.0
+xs = lower_bound:0.1:upper_bound
+x = sample(n_samples, lower_bound, upper_bound, SobolSample())
+y = f.(x)
+gp_surrogate = AbstractGPSurrogate(x,y)
+@show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, gp_surrogate, SobolSample())
+```
+Plotting the function and the sampled points: 
+
+```@example abstractgps_tutorial_optimization
+scatter(gp_surrogate.x, gp_surrogate.y, label="Sampled points", ylims=(-1.0, 5.0), legend=:top)
+plot!(xs, gp_surrogate.(xs), label="Surrogate function", ribbon=p->std_error_at_point(gp_surrogate, p), legend=:top)
+plot!(xs, f.(xs), label="True function", legend=:top)
+```
+
+## ND Example
+
+```@example abstractgps_tutorialnd
+using Plots
+default(c=:matter, legend=false, xlabel="x", ylabel="y")
+using Surrogates 
+
+hypot_func = z -> 3*hypot(z...)+1
+n_samples = 50
+lower_bound = [-1.0, -1.0]
+upper_bound = [1.0, 1.0]
+
+xys = sample(n_samples, lower_bound, upper_bound, SobolSample())
+zs = hypot_func.(xys);
+
+x, y = -2:2, -2:2 
+p1 = surface(x, y, (x1,x2) -> hypot_func((x1,x2))) 
+xs = [xy[1] for xy in xys] 
+ys = [xy[2] for xy in xys] 
+scatter!(xs, ys, zs) 
+p2 = contour(x, y, (x1,x2) -> hypot_func((x1,x2)))
+scatter!(xs, ys)
+plot(p1, p2, title="True function")
+```
+Now let's see how our surrogate performs:
+
+```@example abstractgps_tutorialnd
+gp_surrogate = AbstractGPSurrogate(xys, zs)
+p1 = surface(x, y, (x, y) -> gp_surrogate([x y]))
+scatter!(xs, ys, zs, marker_z=zs)
+p2 = contour(x, y, (x, y) -> gp_surrogate([x y]))
+scatter!(xs, ys, marker_z=zs)
+plot(p1, p2, title="Surrogate")
+```
+
+```@example abstractgps_tutorialnd
+@show gp_surrogate((0.2,0.2))
+```
+
+```@example abstractgps_tutorialnd
+@show hypot_func((0.2,0.2))
+```
+
+And this is our log marginal posterior predictive probability:
+```@example abstractgps_tutorialnd
+@show logpdf_surrogate(gp_surrogate)
+```

src/AbstractGP.jl

@@ -0,0 +1,56 @@
+using AbstractGPs
+
+mutable struct AbstractGPStruct{X, Y, GP, GP_P, S} <: AbstractSurrogate
+    x::X
+    y::Y
+    gp::GP 
+    gp_posterior::GP_P
+    Σy::S
+ end
+
+#constructor
+function AbstractGPSurrogate(x, y; gp = GP(Matern52Kernel()), Σy = 0.1)
+    AbstractGPStruct(x, y, gp, posterior(gp(x, Σy),y), Σy)
+ end
+
+# predictor 
+function (g::AbstractGPStruct)(val)
+    return first(mean(g.gp_posterior([val])))
+end
+
+#add point
+#copies of x and y need to be made because we get 
+#"Error: cannot resize array with shared data " if we push! directly to x and y  
+function add_point!(g::AbstractGPStruct, new_x, new_y)
+    if new_x in g.x
+        println("Adding a sample that already exists, cannot build AbstracgGPSurrogate.")
+        return
+    end
+    x_copy = copy(g.x)
+    push!(x_copy, new_x)
+    y_copy = copy(g.y)
+    push!(y_copy, new_y) 
+    updated_posterior = posterior(g.gp(x_copy, g.Σy), y_copy)
+    g.x, g.y, g.gp_posterior = x_copy, y_copy, updated_posterior
+    nothing
+end
+
+"""
+    Returns diagonal elements of cov(f).
+"""
+function var_at_point(g::AbstractGPStruct, val)
+    return var(g.gp_posterior(val))
+end
+
+"""
+
+"""
+function std_error_at_point(g::AbstractGPStruct, val)
+    return sqrt(first(var(g.gp_posterior([val]))))
+end
+
+
+# Log marginal posterior predictive probability.
+function logpdf_surrogate(g::AbstractGPStruct)
+    return logpdf(g.gp_posterior(g.x), g.y)
+end

src/Surrogates.jl

@@ -13,11 +13,9 @@ include("Lobachevsky.jl")
 include("LinearSurrogate.jl")
 include("InverseDistanceSurrogate.jl")
 include("SecondOrderPolynomialSurrogate.jl")
+include("AbstractGP.jl")
 
 function __init__()
-    @require Stheno="8188c328-b5d6-583d-959b-9690869a5511" begin 
-        include("SthenoKriging.jl")
-    end
 
     @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin
         include("NeuralSurrogate.jl")
@@ -110,7 +108,6 @@ export SVMSurrogate
 export NeuralSurrogate
 export InverseDistanceSurrogate
 export SecondOrderPolynomialSurrogate
-export SthenoKriging
 export Wendland
 export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure
 export LobachevskyStructure, NeuralStructure, RandomForestStructure, SecondOrderPolynomialStructure
@@ -120,4 +117,5 @@ export VariableFidelitySurrogate
 export PolynomialChaosSurrogate
 export EarthSurrogate
 export GEK
+export AbstractGPSurrogate, var_at_point, logpdf_surrogate
 end