## Gradient Enhanced Kriging

Gradient-enhanced Kriging is an extension of kriging which supports gradient information. GEK is usually more accurate than kriging, however, it is not computationally efficient when the number of inputs, the number of sampling points, or both, are high. This is mainly due to the size of the corresponding correlation matrix that increases proportionally with both the number of inputs and the number of sampling points.

Let's have a look to the following function to use Gradient Enhanced Surrogate: $f(x) = sin(x) + 2*x^2$

First of all, we will import Surrogates and Plots packages:

using Surrogates
using Plots
default()

### Sampling

We choose to sample f in 8 points between 0 to 1 using the sample function. The sampling points are chosen using a Sobol sequence, this can be done by passing SobolSample() to the sample function.

n_samples = 10
lower_bound = 2
upper_bound = 10
xs = lower_bound:0.001:upper_bound
x = sample(n_samples, lower_bound, upper_bound, SobolSample())
f(x) = x^3 - 6x^2 + 4x + 12
y1 = f.(x)
der = x -> 3*x^2 - 12*x + 4
y2 = der.(x)
y = vcat(y1,y2)
scatter(x, y1, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top)
plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top)

### Building a surrogate

With our sampled points we can build the Gradient Enhanced Kriging surrogate using the GEK function.

my_gek = GEK(x, y, lower_bound, upper_bound, p = 1.4);
(::GEK{Vector{Float64}, Vector{Float64}, Int64, Int64, Float64, Float64, Float64, Matrix{Float64}, Float64, Matrix{Float64}}) (generic function with 2 methods)
scatter(x, y1, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top)
plot!(f, label="True function",  xlims=(lower_bound, upper_bound), legend=:top)
plot!(my_gek, label="Surrogate function", ribbon=p->std_error_at_point(my_gek, p), xlims=(lower_bound, upper_bound), legend=:top)

## Gradient Enhanced Kriging Surrogate Tutorial (ND)

First of all let's define the function we are going to build a surrogate for.

Now, let's define the function:

function leon(x)
x1 = x
x2 = x
term1 = 100*(x2 - x1^3)^2
term2 = (1 - x1)^2
y = term1 + term2
end
leon (generic function with 1 method)

### Sampling

Let's define our bounds, this time we are working in two dimensions. In particular we want our first dimension x to have bounds 0, 10, and 0, 10 for the second dimension. We are taking 80 samples of the space using Sobol Sequences. We then evaluate our function on all of the sampling points.

n_samples = 80
lower_bound = [0, 0]
upper_bound = [10, 10]
xys = sample(n_samples, lower_bound, upper_bound, SobolSample())
y1 = leon.(xys);
80-element Vector{Float64}:
1577.8677217621007
1.8071765303987595e6
2.1269689354185887e7
19491.21875369508
397142.68190494325
5.1836512480499245e7
7.28543611393035e6
1992.7400718308636
5295.5550515377545
1.2406020272666672e7
⋮
4920.825319683173
5947.404192859795
1.1333865241816055e7
7.268698331698339e7
688180.3567782413
90821.96402000086
3.03575746835075e7
3.53029568501928e6
3643.980667055549

### Building a surrogate

Using the sampled points we build the surrogate, the steps are analogous to the 1-dimensional case.

grad1 = x1 -> 2*(300*(x)^5 - 300*(x)^2*x + x -1)
grad2 = x2 -> 200*(x - (x)^3)
d = 2
n = 10
c = 0
y2 = zeros(eltype(y),n*d)
for i in 1:n
y2[i + c] = grad1(x[i])
y2[i + c + 1] = grad2(x[i])
c = c + 1
end
return y2
end
y = vcat(y1,y2)
100-element Vector{Float64}:
1577.8677217621007
1.8071765303987595e6
2.1269689354185887e7
19491.21875369508
397142.68190494325
5.1836512480499245e7
7.28543611393035e6
1992.7400718308636
5295.5550515377545
1.2406020272666672e7
⋮
200.0
200.0
200.0
200.0
200.0
200.0
200.0
200.0
200.0
my_GEK = GEK(xys, y, lower_bound, upper_bound, p=[1.9, 1.9])
(::GEK{Vector{Tuple{Float64, Float64}}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Float64}, Float64, Matrix{Float64}, Float64, Matrix{Float64}}) (generic function with 2 methods)