Merge branch 'master' of https://github.com/SciML/Surrogates.jl

README.md

@@ -50,6 +50,7 @@ The construction of a surrogate model can be seen as a three-step process:
 - LCBS
 - DYCORS
 - EI
+- SOP
 
 ## Installing Surrogates package
 

docs/make.jl

@@ -5,12 +5,18 @@ makedocs(
     "index.md"
     "Tutorials" => [
         "Basics" => "tutorials.md",
-        "RandomForestSurrogate" => "randomforest.md",
         "Radials" => "radials.md",
         "Kriging" => "kriging.md",
         "Lobachesky" => "lobachesky.md",
-        "LinearSurrogate" => "LinearSurrogate.md",
-        "InverseDistance" => "InverseDistance.md"
+        "Linear" => "LinearSurrogate.md",
+        "InverseDistance" => "InverseDistance.md",
+        "RandomForest" => "randomforest.md",
+        "SecondOrderPolynomial" => "secondorderpoly.md",
+        "NeuralSurrogate" => "neural.md",
+        "Wendland" => "wendland.md",
+        "Polynomial Chaos" => "polychaos.md",
+        "Variable Fidelity" => "variablefidelity.md",
+        "Gradient Enhanced Kriging" => "gek.md"
         ]
     "User guide" => [
         "Samples" => "samples.md",
@@ -25,7 +31,10 @@ makedocs(
         "Cantilever beam" => "cantilever.md",
         "Water Flow function" => "water_flow.md",
         "Welded beam function" => "welded_beam.md",
-        "Branin Function" => "BraninFunction.md"
+        "Branin function" => "BraninFunction.md",
+        "Ackley function" => "ackley.md",
+        "Gramacy & Lee Function" => "gramacylee.md",
+        "Salustowicz Benchmark" => "Salustowicz.md"
         ]
     "Contributing" => "contributing.md"
     ]

docs/src/Salustowicz.md

@@ -0,0 +1,66 @@
+## Salustowicz Benchmark Function
+
+The true underlying function HyGP had to approximate is the 1D Salustowicz function. The function can be evaluated in the given domain:
+``x \in [0, 10]``.
+
+The Salustowicz benchmark function is as follows:
+
+``f(x) = e^(-x) * x^3 * cos(x) * sin(x) * (cos(x) * sin(x)*sin(x) - 1)``
+
+Let's import these two packages  `Surrogates` and `Plots`:
+
+```@example salustowicz1D
+using Surrogates
+using Plots
+default()
+```
+
+Now, let's define our objective function:
+
+```@example salustowicz1D
+function salustowicz(x)
+    term1 = 2.72^(-x) * x^3 * cos(x) * sin(x);
+    term2 = (cos(x) * sin(x)*sin(x) - 1);
+    y = term1 * term2;
+end
+```
+
+Let's sample f in 30 points between 0 and 10 using the `sample` function. The sampling points are chosen using a Sobol Sample, this can be done by passing `SobolSample()` to the `sample` function.
+
+```@example salustowicz1D
+n_samples = 30
+lower_bound = 0
+upper_bound = 10
+num_round = 2
+x = sample(n_samples, lower_bound, upper_bound, SobolSample())
+y = salustowicz.(x)
+xs = lower_bound:0.001:upper_bound
+scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top)
+plot!(xs, salustowicz.(xs), label="True function", legend=:top)
+```
+
+Now, let's fit Salustowicz Function with different Surrogates:
+
+```@example salustowicz1D
+InverseDistance = InverseDistanceSurrogate(x, y, lower_bound, upper_bound)
+randomforest_surrogate = RandomForestSurrogate(x ,y ,lower_bound, upper_bound, num_round = 2)
+lobachevsky_surrogate = LobacheskySurrogate(x, y, lower_bound, upper_bound, alpha = 2.0, n = 6)
+scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:topright)
+plot!(xs, salustowicz.(xs), label="True function", legend=:topright)
+plot!(xs, InverseDistance.(xs), label="InverseDistanceSurrogate", legend=:topright)
+plot!(xs, randomforest_surrogate.(xs), label="RandomForest", legend=:topright)
+plot!(xs, lobachevsky_surrogate.(xs), label="Lobachesky", legend=:topright)
+```
+
+Not's let's see Kriging Surrogate with different hyper parameter:
+
+```@example salustowicz1D
+kriging_surrogate1 = Kriging(x, y, lower_bound, upper_bound, p=0.9);
+kriging_surrogate2 = Kriging(x, y, lower_bound, upper_bound, p=1.5);
+kriging_surrogate3 = Kriging(x, y, lower_bound, upper_bound, p=1.9);
+scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:topright)
+plot!(xs, salustowicz.(xs), label="True function", legend=:topright)
+plot!(xs, kriging_surrogate1.(xs), label="kriging_surrogate1", ribbon=p->std_error_at_point(kriging_surrogate1, p), legend=:topright)
+plot!(xs, kriging_surrogate2.(xs), label="kriging_surrogate2", ribbon=p->std_error_at_point(kriging_surrogate2, p), legend=:topright)
+plot!(xs, kriging_surrogate3.(xs), label="kriging_surrogate3", ribbon=p->std_error_at_point(kriging_surrogate3, p), legend=:topright)
+```

docs/src/ackley.md

@@ -0,0 +1,68 @@
+# Ackley function
+
+The Ackley function is defined as:
+``f(x) = -a*exp(-b\sqrt{\frac{1}{d}\sum_{i=1}^d x_i^2}) - exp(\frac{1}{d} \sum_{i=1}^d cos(cx_i)) + a + exp(1)``
+Usually the recommended values are: ``a =  20``, ``b = 0.2`` and ``c =  2\pi``
+
+Let's see the 1D case.
+
+```@example ackley
+using Surrogates
+using Plots
+default()
+```
+
+Now, let's define the `Ackley` function:
+
+```@example ackley
+function ackley(x)
+    a, b, c = 20.0, -0.2, 2.0*π
+    len_recip = inv(length(x))
+    sum_sqrs = zero(eltype(x))
+    sum_cos = sum_sqrs
+    for i in x
+        sum_cos += cos(c*i)
+        sum_sqrs += i^2
+    end
+    return (-a * exp(b * sqrt(len_recip*sum_sqrs)) -
+            exp(len_recip*sum_cos) + a + 2.71)
+end
+```
+
+
+```@example ackley
+n = 100
+lb = -32.768
+ub = 32.768
+x = sample(n, lb, ub, SobolSample())
+y = ackley.(x)
+xs = lb:0.001:ub
+scatter(x, y, label="Sampled points", xlims=(lb, ub), ylims=(0,30), legend=:top)
+plot!(xs, ackley.(xs), label="True function", legend=:top)
+```
+
+```@example ackley
+my_rad = RadialBasis(x, y, lb, ub)
+my_krig = Kriging(x, y, lb, ub)
+my_loba = LobacheskySurrogate(x, y, lb, ub)
+```
+
+```@example ackley
+scatter(x, y, label="Sampled points", xlims=(lb, ub), ylims=(0, 30), legend=:top)
+plot!(xs, ackley.(xs), label="True function", legend=:top)
+plot!(xs, my_rad.(xs), label="Polynomial expansion", legend=:top)
+plot!(xs, my_krig.(xs), label="Lobachesky", legend=:top)
+plot!(xs, my_loba.(xs), label="Kriging", legend=:top)
+```
+
+The fit looks good. Let's now see if we are able to find the minimum value using
+optimization methods:
+
+```@example ackley
+surrogate_optimize(ackley,DYCORS(),lb,ub,my_rad,UniformSample())
+scatter(x, y, label="Sampled points", xlims=(lb, ub), ylims=(0, 30), legend=:top)
+plot!(xs, ackley.(xs), label="True function", legend=:top)
+plot!(xs, my_rad.(xs), label="Radial basis optimized", legend=:top)
+```
+
+The DYCORS methods successfully finds the minimum.

docs/src/gek.md

@@ -0,0 +1,126 @@
+## 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:
+
+```@example GEK1D
+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.
+
+```@example GEK1D
+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.
+
+```@example GEK1D
+my_gek = GEK(x, y, lower_bound, upper_bound, p = 1.4);
+```
+```@example @GEK1D
+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.
+
+```@example GEK_ND
+using Plots # hide
+default(c=:matter, legend=false, xlabel="x", ylabel="y") # hide
+using Surrogates # hide
+```
+
+Now, let's define the function:
+
+```@example GEK_ND
+function leon(x)
+      x1 = x[1]
+      x2 = x[2]
+      term1 = 100*(x2 - x1^3)^2
+      term2 = (1 - x1)^2
+      y = term1 + term2
+end
+```
+
+### 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.
+
+```@example GEK_ND
+n_samples = 80
+lower_bound = [0, 0]
+upper_bound = [10, 10]
+xys = sample(n_samples, lower_bound, upper_bound, SobolSample())
+y1 = leon.(xys);
+```
+
+```@example GEK_ND
+x, y = 0:10, 0:10 # hide
+p1 = surface(x, y, (x1,x2) -> leon((x1,x2))) # hide
+xs = [xy[1] for xy in xys] # hide
+ys = [xy[2] for xy in xys] # hide
+scatter!(xs, ys, y1) # hide
+p2 = contour(x, y, (x1,x2) -> leon((x1,x2))) # hide
+scatter!(xs, ys) # hide
+plot(p1, p2, title="True function") # hide
+```
+
+### Building a surrogate
+Using the sampled points we build the surrogate, the steps are analogous to the 1-dimensional case.
+
+```@example GEK_ND
+grad1 = x1 -> 2*(300*(x[1])^5 - 300*(x[1])^2*x[2] + x[1] -1)
+grad2 = x2 -> 200*(x[2] - (x[1])^3)
+d = 2
+n = 10
+function create_grads(n, d, grad1, grad2, y)
+      c = 0
+      y2 = zeros(eltype(y[1]),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
+y2 = create_grads(n, d, grad2, grad2, y)
+y = vcat(y1,y2)
+```
+
+```@example GEK_ND
+my_GEK = GEK(xys, y, lower_bound, upper_bound, p=[1.9, 1.9])
+```
+
+```@example GEK_ND
+p1 = surface(x, y, (x, y) -> my_GEK([x y])) # hide
+scatter!(xs, ys, y1, marker_z=y1) # hide
+p2 = contour(x, y, (x, y) -> my_GEK([x y])) # hide
+scatter!(xs, ys, marker_z=y1) # hide
+plot(p1, p2, title="Surrogate") # hide
+```

docs/src/gramacylee.md

@@ -0,0 +1,51 @@
+## Gramacy & Lee Function
+
+Gramacy & Lee Function is a continues function. It is not convex. The function is defined on 1-dimensional space. It is an unimodal. The function can be defined on any input domain but it is usually evaluated on
+``x \in [-0.5, 2.5]``.
+
+The Gramacy & Lee is as follows:
+``f(x) = \frac{sin(10\pi x)}{2x} + (x-1)^4``.
+
+Let's import these two packages `Surrogates` and `Plots`:
+
+```@example gramacylee1D
+using Surrogates
+using Plots
+default()
+```
+
+Now, let's define our objective function:
+
+```@example gramacylee1D
+function gramacylee(x)
+    term1 = sin(10*pi*x) / 2*x;
+    term2 = (x - 1)^4;
+    y = term1 + term2;
+end
+```
+
+Let's sample f in 25 points between -0.5 and 2.5 using the `sample` function. The sampling points are chosen using a Sobol Sample, this can be done by passing `SobolSample()` to the `sample` function.
+
+```@example gramacylee1D
+n = 25
+lower_bound = -0.5
+upper_bound = 2.5
+x = sample(n, lower_bound, upper_bound, SobolSample())
+y = gramacylee.(x)
+xs = lower_bound:0.001:upper_bound
+scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound), ylims=(-5, 20), legend=:top)
+plot!(xs, gramacylee.(xs), label="True function", legend=:top)
+```
+
+Now, let's fit Gramacy & Lee Function with different Surrogates:
+
+```@example gramacylee1D
+my_pol = PolynomialChaosSurrogate(x, y, lower_bound, upper_bound)
+loba_1 = LobacheskySurrogate(x, y, lower_bound, upper_bound)
+krig = Kriging(x, y, lower_bound, upper_bound)
+scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound), ylims=(-5, 20), legend=:top)
+plot!(xs, gramacylee.(xs), label="True function", legend=:top)
+plot!(xs, my_pol.(xs), label="Polynomial expansion", legend=:top)
+plot!(xs, loba_1.(xs), label="Lobachesky", legend=:top)
+plot!(xs, krig.(xs), label="Kriging", legend=:top)
+```