From 7d22d8ebc26bbe4c5b139048f48202696ba9152d Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Sat, 24 Feb 2024 01:45:09 +0100 Subject: [PATCH] reapply formatter --- .JuliaFormatter.toml | 4 +- CONTRIBUTING.md | 28 +- LICENSE.md | 6 +- README.md | 53 +-- docs/pages.jl | 78 ++--- docs/src/BraninFunction.md | 188 +++++----- docs/src/ImprovedBraninFunction.md | 18 +- docs/src/InverseDistance.md | 304 ++++++++-------- docs/src/LinearSurrogate.md | 299 ++++++++-------- docs/src/Salustowicz.md | 134 ++++---- docs/src/abstractgps.md | 69 ++-- docs/src/ackley.md | 32 +- docs/src/cantilever.md | 45 +-- docs/src/gek.md | 251 +++++++------- docs/src/gekpls.md | 83 ++--- docs/src/gramacylee.md | 22 +- docs/src/index.md | 80 +++-- docs/src/kriging.md | 80 +++-- docs/src/lobachevsky.md | 74 ++-- docs/src/lp.md | 36 +- docs/src/moe.md | 26 +- docs/src/multi_objective_opt.md | 38 +- docs/src/neural.md | 36 +- docs/src/optimizations.md | 18 +- docs/src/parallel.md | 49 +-- docs/src/polychaos.md | 27 +- docs/src/radials.md | 325 +++++++++--------- docs/src/randomforest.md | 76 ++-- docs/src/rosenbrock.md | 42 ++- docs/src/samples.md | 25 +- docs/src/secondorderpoly.md | 22 +- docs/src/sphere_function.md | 49 +-- docs/src/surrogate.md | 38 +- docs/src/tensor_prod.md | 21 +- docs/src/tutorials.md | 49 +-- docs/src/variablefidelity.md | 19 +- docs/src/water_flow.md | 24 +- docs/src/welded_beam.md | 26 +- docs/src/wendland.md | 11 +- lib/SurrogatesAbstractGPs/LICENSE.md | 6 +- lib/SurrogatesFlux/LICENSE.md | 6 +- lib/SurrogatesFlux/src/SurrogatesFlux.jl | 7 +- lib/SurrogatesMOE/LICENSE.md | 6 +- lib/SurrogatesMOE/src/SurrogatesMOE.jl | 31 +- lib/SurrogatesMOE/test/runtests.jl | 14 +- lib/SurrogatesPolyChaos/LICENSE.md | 6 +- lib/SurrogatesRandomForest/LICENSE.md | 6 +- .../src/SurrogatesRandomForest.jl | 1 - lib/SurrogatesSVM/LICENSE.md | 6 +- src/Earth.jl | 6 +- src/GEK.jl | 2 +- src/GEKPLS.jl | 146 ++++---- src/InverseDistanceSurrogate.jl | 5 +- src/Kriging.jl | 44 +-- src/LinearSurrogate.jl | 1 - src/Lobachevsky.jl | 9 +- src/Optimization.jl | 16 +- src/Radials.jl | 13 +- src/Sampling.jl | 5 +- src/Surrogates.jl | 12 +- src/VariableFidelity.jl | 2 - 61 files changed, 1666 insertions(+), 1489 deletions(-) diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 453925c3f..959ad88a6 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -1 +1,3 @@ -style = "sciml" \ No newline at end of file +style = "sciml" +format_markdown = true +format_docstrings = true \ No newline at end of file diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 6d1b91fa2..c4d2791f0 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,21 +1,13 @@ # Contributing -1. Fork the repository on github. (Click the `Fork` button in the top-right corner) - -2. Clone the repository you have just forked. `git clone https://github.com/YOUR_USERNAME/Surrogates.jl.git` - -3. `cd Surrogates.jl` Enter the repository's directory. - -4. `julia` Open the Julia REPL. - -5. `]` Enter package mode - -6. Activate the local environment `activate .` - -7. Install the dependencies. `instantiate` - -8. Perform your edits (Atom with Juno, or VSCode with the Julia plugin are good editor choices) - -9. Stage, Commit, and Push your changes - + 1. Fork the repository on github. (Click the `Fork` button in the top-right corner) + + 2. Clone the repository you have just forked. `git clone https://github.com/YOUR_USERNAME/Surrogates.jl.git` + 3. `cd Surrogates.jl` Enter the repository's directory. + 4. `julia` Open the Julia REPL. + 5. `]` Enter package mode + 6. Activate the local environment `activate .` + 7. Install the dependencies. `instantiate` + 8. Perform your edits (Atom with Juno, or VSCode with the Julia plugin are good editor choices) + 9. Stage, Commit, and Push your changes 10. [Open a Pull Request](https://help.github.com/en/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request-from-a-fork) diff --git a/LICENSE.md b/LICENSE.md index 62199eff7..27dd257c0 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,7 +1,7 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > Copyright (c) 2019-20: Ludovico Bessi, Julia Computing. -> +> > Permission is hereby granted, free of charge, to any person obtaining > a copy of this software and associated documentation files (the > "Software"), to deal in the Software without restriction, including @@ -9,10 +9,10 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > distribute, sublicense, and/or sell copies of the Software, and to > permit persons to whom the Software is furnished to do so, subject to > the following conditions: -> +> > The above copyright notice and this permission notice shall be > included in all copies or substantial portions of the Software. -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. diff --git a/README.md b/README.md index 40a0427fb..6e065a6df 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ [![codecov](https://codecov.io/gh/SciML/Surrogates.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/SciML/Surrogates.jl) [![Build Status](https://github.com/SciML/Surrogates.jl/workflows/CI/badge.svg)](https://github.com/SciML/Surrogates.jl/actions?query=workflow%3ACI) -[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac) +[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac) [![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle) A surrogate model is an approximation method that mimics the behavior of a computationally @@ -14,39 +14,40 @@ expensive simulation. In more mathematical terms: suppose we are attempting to o `f(p)`, but each calculation of `f` is very expensive. It may be the case we need to solve a PDE for each point or use advanced numerical linear algebra machinery, which is usually costly. The idea is then to develop a surrogate model `g` which approximates `f` by training on previous data collected from evaluations of `f`. The construction of a surrogate model can be seen as a three-step process: -1. Sample selection -2. Construction of the surrogate model -3. Surrogate optimization + 1. Sample selection + 2. Construction of the surrogate model + 3. Surrogate optimization Sampling can be done through [QuasiMonteCarlo.jl](https://github.com/SciML/QuasiMonteCarlo.jl), all the functions available there can be used in Surrogates.jl. ## ALL the currently available surrogate models: -- Kriging -- Kriging using Stheno -- Radial Basis -- Wendland -- Linear -- Second Order Polynomial -- Support Vector Machines (Wait for LIBSVM resolution) -- Neural Networks -- Random Forests -- Lobachevsky -- Inverse-distance -- Polynomial expansions -- Variable fidelity -- Mixture of experts (Waiting GaussianMixtures package to work on v1.5) -- Earth -- Gradient Enhanced Kriging + - Kriging + - Kriging using Stheno + - Radial Basis + - Wendland + - Linear + - Second Order Polynomial + - Support Vector Machines (Wait for LIBSVM resolution) + - Neural Networks + - Random Forests + - Lobachevsky + - Inverse-distance + - Polynomial expansions + - Variable fidelity + - Mixture of experts (Waiting GaussianMixtures package to work on v1.5) + - Earth + - Gradient Enhanced Kriging ## ALL the currently available optimization methods: -- SRBF -- LCBS -- DYCORS -- EI -- SOP -- Multi-optimization: SMB and RTEA + - SRBF + - LCBS + - DYCORS + - EI + - SOP + - Multi-optimization: SMB and RTEA + ## Installing Surrogates package ```julia diff --git a/docs/pages.jl b/docs/pages.jl index 24d4d3a79..b40c7df1b 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,40 +1,40 @@ pages = ["index.md" - "Tutorials" => [ - "Basics" => "tutorials.md", - "Radials" => "radials.md", - "Kriging" => "kriging.md", - "Gaussian Process" => "abstractgps.md", - "Lobachevsky" => "lobachevsky.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", - "GEKPLS" => "gekpls.md", - "MOE" => "moe.md", - "Parallel Optimization" => "parallel.md", -] - "User guide" => [ - "Samples" => "samples.md", - "Surrogates" => "surrogate.md", - "Optimization" => "optimizations.md", -] - "Benchmarks" => [ - "Sphere function" => "sphere_function.md", - "Lp norm" => "lp.md", - "Rosenbrock" => "rosenbrock.md", - "Tensor product" => "tensor_prod.md", - "Cantilever beam" => "cantilever.md", - "Water Flow function" => "water_flow.md", - "Welded beam function" => "welded_beam.md", - "Branin function" => "BraninFunction.md", - "Improved Branin function" => "ImprovedBraninFunction.md", - "Ackley function" => "ackley.md", - "Gramacy & Lee Function" => "gramacylee.md", - "Salustowicz Benchmark" => "Salustowicz.md", - "Multi objective optimization" => "multi_objective_opt.md", -]] + "Tutorials" => [ + "Basics" => "tutorials.md", + "Radials" => "radials.md", + "Kriging" => "kriging.md", + "Gaussian Process" => "abstractgps.md", + "Lobachevsky" => "lobachevsky.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", + "GEKPLS" => "gekpls.md", + "MOE" => "moe.md", + "Parallel Optimization" => "parallel.md" + ] + "User guide" => [ + "Samples" => "samples.md", + "Surrogates" => "surrogate.md", + "Optimization" => "optimizations.md" + ] + "Benchmarks" => [ + "Sphere function" => "sphere_function.md", + "Lp norm" => "lp.md", + "Rosenbrock" => "rosenbrock.md", + "Tensor product" => "tensor_prod.md", + "Cantilever beam" => "cantilever.md", + "Water Flow function" => "water_flow.md", + "Welded beam function" => "welded_beam.md", + "Branin function" => "BraninFunction.md", + "Improved Branin function" => "ImprovedBraninFunction.md", + "Ackley function" => "ackley.md", + "Gramacy & Lee Function" => "gramacylee.md", + "Salustowicz Benchmark" => "Salustowicz.md", + "Multi objective optimization" => "multi_objective_opt.md" + ]] diff --git a/docs/src/BraninFunction.md b/docs/src/BraninFunction.md index c58abe12a..49e98e0fc 100644 --- a/docs/src/BraninFunction.md +++ b/docs/src/BraninFunction.md @@ -1,92 +1,96 @@ -# Branin Function - -The Branin function is commonly used as a test function for metamodelling in computer experiments, especially in the context of optimization. - -The expression of the Branin Function is given as: -``f(x) = (x_2 - \frac{5.1}{4\pi^2}x_1^{2} + \frac{5}{\pi}x_1 - 6)^2 + 10(1-\frac{1}{8\pi})\cos(x_1) + 10`` - -where ``x = (x_1, x_2)`` with ``-5\leq x_1 \leq 10, 0 \leq x_2 \leq 15`` - -First of all, we will import these two packages: `Surrogates` and `Plots`. - -```@example BraninFunction -using Surrogates -using Plots -default() -``` - -Now, let's define our objective function: - -```@example BraninFunction -function branin(x) - x1 = x[1] - x2 = x[2] - b = 5.1 / (4*pi^2); - c = 5/pi; - r = 6; - a = 1; - s = 10; - t = 1 / (8*pi); - term1 = a * (x2 - b*x1^2 + c*x1 - r)^2; - term2 = s*(1-t)*cos(x1); - y = term1 + term2 + s; -end -``` -Now, let's plot it: - -```@example BraninFunction -n_samples = 80 -lower_bound = [-5, 0] -upper_bound = [10,15] -xys = sample(n_samples, lower_bound, upper_bound, SobolSample()) -zs = branin.(xys); -x, y = -5.00:10.00, 0.00:15.00 -p1 = surface(x, y, (x1,x2) -> branin((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) -> branin((x1,x2))) -scatter!(xs, ys) -plot(p1, p2, title="True function") -``` - -Now it's time to try fitting different surrogates, and then we will plot them. -We will have a look at the radial basis surrogate `Radial Basis Surrogate`. : - -```@example BraninFunction -radial_surrogate = RadialBasis(xys, zs, lower_bound, upper_bound) -``` - -```@example BraninFunction -p1 = surface(x, y, (x, y) -> radial_surrogate([x y])) -scatter!(xs, ys, zs, marker_z=zs) -p2 = contour(x, y, (x, y) -> radial_surrogate([x y])) -scatter!(xs, ys, marker_z=zs) -plot(p1, p2, title="Radial Surrogate") -``` - -Now, we will have a look at `Inverse Distance Surrogate`: -```@example BraninFunction -InverseDistance = InverseDistanceSurrogate(xys, zs, lower_bound, upper_bound) -``` - -```@example BraninFunction -p1 = surface(x, y, (x, y) -> InverseDistance([x y])) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide -p2 = contour(x, y, (x, y) -> InverseDistance([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2, title="Inverse Distance Surrogate") # hide -``` - -Now, let's talk about `Lobachevsky Surrogate`: -```@example BraninFunction -Lobachevsky = LobachevskySurrogate(xys, zs, lower_bound, upper_bound, alpha = [2.8,2.8], n=8) -``` - -```@example BraninFunction -p1 = surface(x, y, (x, y) -> Lobachevsky([x y])) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide -p2 = contour(x, y, (x, y) -> Lobachevsky([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2, title="Lobachevsky Surrogate") # hide -``` +# Branin Function + +The Branin function is commonly used as a test function for metamodelling in computer experiments, especially in the context of optimization. + +The expression of the Branin Function is given as: +``f(x) = (x_2 - \frac{5.1}{4\pi^2}x_1^{2} + \frac{5}{\pi}x_1 - 6)^2 + 10(1-\frac{1}{8\pi})\cos(x_1) + 10`` + +where ``x = (x_1, x_2)`` with ``-5\leq x_1 \leq 10, 0 \leq x_2 \leq 15`` + +First of all, we will import these two packages: `Surrogates` and `Plots`. + +```@example BraninFunction +using Surrogates +using Plots +default() +``` + +Now, let's define our objective function: + +```@example BraninFunction +function branin(x) + x1 = x[1] + x2 = x[2] + b = 5.1 / (4 * pi^2) + c = 5 / pi + r = 6 + a = 1 + s = 10 + t = 1 / (8 * pi) + term1 = a * (x2 - b * x1^2 + c * x1 - r)^2 + term2 = s * (1 - t) * cos(x1) + y = term1 + term2 + s +end +``` + +Now, let's plot it: + +```@example BraninFunction +n_samples = 80 +lower_bound = [-5, 0] +upper_bound = [10, 15] +xys = sample(n_samples, lower_bound, upper_bound, SobolSample()) +zs = branin.(xys); +x, y = -5.00:10.00, 0.00:15.00 +p1 = surface(x, y, (x1, x2) -> branin((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) -> branin((x1, x2))) +scatter!(xs, ys) +plot(p1, p2, title = "True function") +``` + +Now it's time to try fitting different surrogates, and then we will plot them. +We will have a look at the radial basis surrogate `Radial Basis Surrogate`. : + +```@example BraninFunction +radial_surrogate = RadialBasis(xys, zs, lower_bound, upper_bound) +``` + +```@example BraninFunction +p1 = surface(x, y, (x, y) -> radial_surrogate([x y])) +scatter!(xs, ys, zs, marker_z = zs) +p2 = contour(x, y, (x, y) -> radial_surrogate([x y])) +scatter!(xs, ys, marker_z = zs) +plot(p1, p2, title = "Radial Surrogate") +``` + +Now, we will have a look at `Inverse Distance Surrogate`: + +```@example BraninFunction +InverseDistance = InverseDistanceSurrogate(xys, zs, lower_bound, upper_bound) +``` + +```@example BraninFunction +p1 = surface(x, y, (x, y) -> InverseDistance([x y])) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide +p2 = contour(x, y, (x, y) -> InverseDistance([x y])) # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2, title = "Inverse Distance Surrogate") # hide +``` + +Now, let's talk about `Lobachevsky Surrogate`: + +```@example BraninFunction +Lobachevsky = LobachevskySurrogate( + xys, zs, lower_bound, upper_bound, alpha = [2.8, 2.8], n = 8) +``` + +```@example BraninFunction +p1 = surface(x, y, (x, y) -> Lobachevsky([x y])) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide +p2 = contour(x, y, (x, y) -> Lobachevsky([x y])) # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2, title = "Lobachevsky Surrogate") # hide +``` diff --git a/docs/src/ImprovedBraninFunction.md b/docs/src/ImprovedBraninFunction.md index 002a5ef2f..1a65aede5 100644 --- a/docs/src/ImprovedBraninFunction.md +++ b/docs/src/ImprovedBraninFunction.md @@ -10,17 +10,17 @@ To enhance the Branin function, changes were made to introduce irregularities, v function improved_branin(x, time_step) x1 = x[1] x2 = x[2] - b = 5.1 / (4*pi^2) - c = 5/pi + b = 5.1 / (4 * pi^2) + c = 5 / pi r = 6 a = 1 s = 10 - t = 1 / (8*pi) - + t = 1 / (8 * pi) + # Adding noise to the function's output noise = randn() * time_step # Simulating time-varying noise - term1 = a * (x2 - b*x1^2 + c*x1 - r)^2 - term2 = s*(1-t)*cos(x1 + noise) # Introducing dynamic component + term1 = a * (x2 - b * x1^2 + c * x1 - r)^2 + term2 = s * (1 - t) * cos(x1 + noise) # Introducing dynamic component y = term1 + term2 + s end ``` @@ -44,8 +44,8 @@ x, y = -5.00:10.00, 0.00:15.00 xs = [xy[1] for xy in xys] ys = [xy[2] for xy in xys] p1 = surface(x, y, (x, y) -> radial_surrogate([x, y])) -scatter!(xs, ys, marker_z=zs) +scatter!(xs, ys, marker_z = zs) p2 = contour(x, y, (x, y) -> radial_surrogate([x, y])) -scatter!(xs, ys, marker_z=zs) -plot(p1, p2, title="Radial Surrogate") +scatter!(xs, ys, marker_z = zs) +plot(p1, p2, title = "Radial Surrogate") ``` diff --git a/docs/src/InverseDistance.md b/docs/src/InverseDistance.md index 2ed3115d1..814b60c1b 100644 --- a/docs/src/InverseDistance.md +++ b/docs/src/InverseDistance.md @@ -1,151 +1,153 @@ -The **Inverse Distance Surrogate** is an interpolating method, and in this method, the unknown points are calculated with a weighted average of the sampling points. This model uses the inverse distance between the unknown and training points to predict the unknown point. We do not need to fit this model because the response of an unknown point x is computed with respect to the distance between x and the training points. - -Let's optimize the following function to use Inverse Distance Surrogate: - -$f(x) = sin(x) + sin(x)^2 + sin(x)^3$. - -First of all, we have to import these two packages: `Surrogates` and `Plots`. - -```@example Inverse_Distance1D -using Surrogates -using Plots -default() -``` - - -### Sampling - -We choose to sample f in 25 points between 0 and 10 using the `sample` function. The sampling points are chosen using a Low Discrepancy, this can be done by passing `HaltonSample()` to the `sample` function. - -```@example Inverse_Distance1D -f(x) = sin(x) + sin(x)^2 + sin(x)^3 - -n_samples = 25 -lower_bound = 0.0 -upper_bound = 10.0 -x = sample(n_samples, lower_bound, upper_bound, HaltonSample()) -y = f.(x) - -scatter(x, y, 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 - -```@example Inverse_Distance1D -InverseDistance = InverseDistanceSurrogate(x, y, lower_bound, upper_bound) -add_point!(InverseDistance, 5.0, f(5.0)) -add_point!(InverseDistance, [5.1,5.2], [f(5.1),f(5.2)]) -prediction = InverseDistance(5.0) -``` - -Now, we will simply plot `InverseDistance`: - -```@example Inverse_Distance1D -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) -plot!(InverseDistance, label="Surrogate function", xlims=(lower_bound, upper_bound), legend=:top) -``` - - -## Optimizing - -Having built a surrogate, we can now use it to search for minima in our original function `f`. - -To optimize using our surrogate we call `surrogate_optimize` method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique. - -```@example Inverse_Distance1D -@show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, InverseDistance, SobolSample()) -scatter(x, y, label="Sampled points", legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) -plot!(InverseDistance, label="Surrogate function", xlims=(lower_bound, upper_bound), legend=:top) -``` - - -## Inverse Distance Surrogate Tutorial (ND): - -First of all we will define the `Schaffer` function we are going to build a surrogate for. Notice, how its argument is a vector of numbers, one for each coordinate, and its output is a scalar. - -```@example Inverse_DistanceND -using Plots # hide -default(c=:matter, legend=false, xlabel="x", ylabel="y") # hide -using Surrogates # hide - -function schaffer(x) - x1=x[1] - x2=x[2] - fact1 = (sin(x1^2-x2^2))^2 - 0.5; - fact2 = (1 + 0.001*(x1^2+x2^2))^2; - y = 0.5 + fact1/fact2; -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 `-5, 10`, and `0, 15` for the second dimension. We are taking 60 samples of the space using Sobol Sequences. We then evaluate our function on all the sampling points. - -```@example Inverse_DistanceND -n_samples = 60 -lower_bound = [-5.0, 0.0] -upper_bound = [10.0, 15.0] - -xys = sample(n_samples, lower_bound, upper_bound, SobolSample()) -zs = schaffer.(xys); -``` - -```@example Inverse_DistanceND -x, y = -5:10, 0:15 # hide -p1 = surface(x, y, (x1,x2) -> schaffer((x1,x2))) # hide -xs = [xy[1] for xy in xys] # hide -ys = [xy[2] for xy in xys] # hide -scatter!(xs, ys, zs) # hide -p2 = contour(x, y, (x1,x2) -> schaffer((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 Inverse_DistanceND -InverseDistance = InverseDistanceSurrogate(xys, zs, lower_bound, upper_bound) -``` - -```@example Inverse_DistanceND -p1 = surface(x, y, (x, y) -> InverseDistance([x y])) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide -p2 = contour(x, y, (x, y) -> InverseDistance([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2, title="Surrogate") # hide -``` - - -### Optimizing -With our surrogate, we can now search for the minima of the function. - -Notice how the new sampled points, which were created during the optimization process, are appended to the `xys` array. -This is why its size changes. - -```@example Inverse_DistanceND -size(xys) -``` -```@example Inverse_DistanceND -surrogate_optimize(schaffer, SRBF(), lower_bound, upper_bound, InverseDistance, SobolSample(), maxiters=10) -``` -```@example Inverse_DistanceND -size(xys) -``` - -```@example Inverse_DistanceND -p1 = surface(x, y, (x, y) -> InverseDistance([x y])) # hide -xs = [xy[1] for xy in xys] # hide -ys = [xy[2] for xy in xys] # hide -zs = schaffer.(xys) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide -p2 = contour(x, y, (x, y) -> InverseDistance([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2) # hide -``` +The **Inverse Distance Surrogate** is an interpolating method, and in this method, the unknown points are calculated with a weighted average of the sampling points. This model uses the inverse distance between the unknown and training points to predict the unknown point. We do not need to fit this model because the response of an unknown point x is computed with respect to the distance between x and the training points. + +Let's optimize the following function to use Inverse Distance Surrogate: + +$f(x) = sin(x) + sin(x)^2 + sin(x)^3$. + +First of all, we have to import these two packages: `Surrogates` and `Plots`. + +```@example Inverse_Distance1D +using Surrogates +using Plots +default() +``` + +### Sampling + +We choose to sample f in 25 points between 0 and 10 using the `sample` function. The sampling points are chosen using a Low Discrepancy, this can be done by passing `HaltonSample()` to the `sample` function. + +```@example Inverse_Distance1D +f(x) = sin(x) + sin(x)^2 + sin(x)^3 + +n_samples = 25 +lower_bound = 0.0 +upper_bound = 10.0 +x = sample(n_samples, lower_bound, upper_bound, HaltonSample()) +y = f.(x) + +scatter(x, y, 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 + +```@example Inverse_Distance1D +InverseDistance = InverseDistanceSurrogate(x, y, lower_bound, upper_bound) +add_point!(InverseDistance, 5.0, f(5.0)) +add_point!(InverseDistance, [5.1, 5.2], [f(5.1), f(5.2)]) +prediction = InverseDistance(5.0) +``` + +Now, we will simply plot `InverseDistance`: + +```@example Inverse_Distance1D +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lower_bound, upper_bound), legend = :top) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) +plot!(InverseDistance, label = "Surrogate function", + xlims = (lower_bound, upper_bound), legend = :top) +``` + +## Optimizing + +Having built a surrogate, we can now use it to search for minima in our original function `f`. + +To optimize using our surrogate we call `surrogate_optimize` method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique. + +```@example Inverse_Distance1D +@show surrogate_optimize( + f, SRBF(), lower_bound, upper_bound, InverseDistance, SobolSample()) +scatter(x, y, label = "Sampled points", legend = :top) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) +plot!(InverseDistance, label = "Surrogate function", + xlims = (lower_bound, upper_bound), legend = :top) +``` + +## Inverse Distance Surrogate Tutorial (ND): + +First of all we will define the `Schaffer` function we are going to build a surrogate for. Notice, how its argument is a vector of numbers, one for each coordinate, and its output is a scalar. + +```@example Inverse_DistanceND +using Plots # hide +default(c = :matter, legend = false, xlabel = "x", ylabel = "y") # hide +using Surrogates # hide + +function schaffer(x) + x1 = x[1] + x2 = x[2] + fact1 = (sin(x1^2 - x2^2))^2 - 0.5 + fact2 = (1 + 0.001 * (x1^2 + x2^2))^2 + y = 0.5 + fact1 / fact2 +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 `-5, 10`, and `0, 15` for the second dimension. We are taking 60 samples of the space using Sobol Sequences. We then evaluate our function on all the sampling points. + +```@example Inverse_DistanceND +n_samples = 60 +lower_bound = [-5.0, 0.0] +upper_bound = [10.0, 15.0] + +xys = sample(n_samples, lower_bound, upper_bound, SobolSample()) +zs = schaffer.(xys); +``` + +```@example Inverse_DistanceND +x, y = -5:10, 0:15 # hide +p1 = surface(x, y, (x1, x2) -> schaffer((x1, x2))) # hide +xs = [xy[1] for xy in xys] # hide +ys = [xy[2] for xy in xys] # hide +scatter!(xs, ys, zs) # hide +p2 = contour(x, y, (x1, x2) -> schaffer((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 Inverse_DistanceND +InverseDistance = InverseDistanceSurrogate(xys, zs, lower_bound, upper_bound) +``` + +```@example Inverse_DistanceND +p1 = surface(x, y, (x, y) -> InverseDistance([x y])) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide +p2 = contour(x, y, (x, y) -> InverseDistance([x y])) # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2, title = "Surrogate") # hide +``` + +### Optimizing + +With our surrogate, we can now search for the minima of the function. + +Notice how the new sampled points, which were created during the optimization process, are appended to the `xys` array. +This is why its size changes. + +```@example Inverse_DistanceND +size(xys) +``` + +```@example Inverse_DistanceND +surrogate_optimize(schaffer, SRBF(), lower_bound, upper_bound, + InverseDistance, SobolSample(), maxiters = 10) +``` + +```@example Inverse_DistanceND +size(xys) +``` + +```@example Inverse_DistanceND +p1 = surface(x, y, (x, y) -> InverseDistance([x y])) # hide +xs = [xy[1] for xy in xys] # hide +ys = [xy[2] for xy in xys] # hide +zs = schaffer.(xys) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide +p2 = contour(x, y, (x, y) -> InverseDistance([x y])) # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2) # hide +``` diff --git a/docs/src/LinearSurrogate.md b/docs/src/LinearSurrogate.md index 2e1ca9de3..d4121f097 100644 --- a/docs/src/LinearSurrogate.md +++ b/docs/src/LinearSurrogate.md @@ -1,146 +1,153 @@ -## Linear Surrogate -Linear Surrogate is a linear approach to modeling the relationship between a scalar response or dependent variable and one or more explanatory variables. We will use Linear Surrogate to optimize following function: - -$f(x) = \sin(x) + \log(x)$ - -First of all we have to import these two packages: `Surrogates` and `Plots`. - -```@example linear_surrogate1D -using Surrogates -using Plots -default() -``` - -### Sampling - -We choose to sample f in 20 points between 0 and 10 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 linear_surrogate1D -f(x) = sin(x) + log(x) -n_samples = 20 -lower_bound = 5.2 -upper_bound = 12.5 -x = sample(n_samples, lower_bound, upper_bound, SobolSample()) -y = f.(x) -scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound)) -plot!(f, label="True function", xlims=(lower_bound, upper_bound)) -``` - -## Building a Surrogate - -With our sampled points, we can build the **Linear Surrogate** using the `LinearSurrogate` function. - -We can simply calculate `linear_surrogate` for any value. - -```@example linear_surrogate1D -my_linear_surr_1D = LinearSurrogate(x, y, lower_bound, upper_bound) -add_point!(my_linear_surr_1D,4.0,7.2) -add_point!(my_linear_surr_1D,[5.0,6.0],[8.3,9.7]) -val = my_linear_surr_1D(5.0) -``` - -Now, we will simply plot `linear_surrogate`: - -```@example linear_surrogate1D -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound)) -plot!(f, label="True function", xlims=(lower_bound, upper_bound)) -plot!(my_linear_surr_1D, label="Surrogate function", xlims=(lower_bound, upper_bound)) -``` - -## Optimizing - -Having built a surrogate, we can now use it to search for minima in our original function `f`. - -To optimize using our surrogate we call `surrogate_optimize` method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique. - -```@example linear_surrogate1D -@show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, my_linear_surr_1D, SobolSample()) -scatter(x, y, label="Sampled points") -plot!(f, label="True function", xlims=(lower_bound, upper_bound)) -plot!(my_linear_surr_1D, label="Surrogate function", xlims=(lower_bound, upper_bound)) -``` - - -## Linear Surrogate tutorial (ND) - -First of all we will define the `Egg Holder` function we are going to build a surrogate for. Notice, one how its argument is a vector of numbers, one for each coordinate, and its output is a scalar. - -```@example linear_surrogateND -using Plots # hide -default(c=:matter, legend=false, xlabel="x", ylabel="y") # hide -using Surrogates # hide - -function egg(x) - x1=x[1] - x2=x[2] - term1 = -(x2+47) * sin(sqrt(abs(x2+x1/2+47))); - term2 = -x1 * sin(sqrt(abs(x1-(x2+47)))); - 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 `-10, 5`, and `0, 15` for the second dimension. We are taking 50 samples of the space using Sobol Sequences. We then evaluate our function on all of the sampling points. - -```@example linear_surrogateND -n_samples = 50 -lower_bound = [-10.0, 0.0] -upper_bound = [5.0, 15.0] - -xys = sample(n_samples, lower_bound, upper_bound, SobolSample()) -zs = egg.(xys); -``` - -```@example linear_surrogateND -x, y = -10:5, 0:15 # hide -p1 = surface(x, y, (x1,x2) -> egg((x1,x2))) # hide -xs = [xy[1] for xy in xys] # hide -ys = [xy[2] for xy in xys] # hide -scatter!(xs, ys, zs) # hide -p2 = contour(x, y, (x1,x2) -> egg((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 linear_surrogateND -my_linear_ND = LinearSurrogate(xys, zs, lower_bound, upper_bound) -``` - -```@example linear_surrogateND -p1 = surface(x, y, (x, y) -> my_linear_ND([x y])) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide -p2 = contour(x, y, (x, y) -> my_linear_ND([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2, title="Surrogate") # hide -``` - -### Optimizing -With our surrogate, we can now search for the minima of the function. - -Notice how the new sampled points, which were created during the optimization process, are appended to the `xys` array. -This is why its size changes. - -```@example linear_surrogateND -size(xys) -``` -```@example linear_surrogateND -surrogate_optimize(egg, SRBF(), lower_bound, upper_bound, my_linear_ND, SobolSample(), maxiters=10) -``` -```@example linear_surrogateND -size(xys) -``` - -```@example linear_surrogateND -p1 = surface(x, y, (x, y) -> my_linear_ND([x y])) # hide -xs = [xy[1] for xy in xys] # hide -ys = [xy[2] for xy in xys] # hide -zs = egg.(xys) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide -p2 = contour(x, y, (x, y) -> my_linear_ND([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2) # hide -``` +## Linear Surrogate + +Linear Surrogate is a linear approach to modeling the relationship between a scalar response or dependent variable and one or more explanatory variables. We will use Linear Surrogate to optimize following function: + +$f(x) = \sin(x) + \log(x)$ + +First of all we have to import these two packages: `Surrogates` and `Plots`. + +```@example linear_surrogate1D +using Surrogates +using Plots +default() +``` + +### Sampling + +We choose to sample f in 20 points between 0 and 10 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 linear_surrogate1D +f(x) = sin(x) + log(x) +n_samples = 20 +lower_bound = 5.2 +upper_bound = 12.5 +x = sample(n_samples, lower_bound, upper_bound, SobolSample()) +y = f.(x) +scatter(x, y, label = "Sampled points", xlims = (lower_bound, upper_bound)) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound)) +``` + +## Building a Surrogate + +With our sampled points, we can build the **Linear Surrogate** using the `LinearSurrogate` function. + +We can simply calculate `linear_surrogate` for any value. + +```@example linear_surrogate1D +my_linear_surr_1D = LinearSurrogate(x, y, lower_bound, upper_bound) +add_point!(my_linear_surr_1D, 4.0, 7.2) +add_point!(my_linear_surr_1D, [5.0, 6.0], [8.3, 9.7]) +val = my_linear_surr_1D(5.0) +``` + +Now, we will simply plot `linear_surrogate`: + +```@example linear_surrogate1D +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lower_bound, upper_bound)) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound)) +plot!(my_linear_surr_1D, label = "Surrogate function", xlims = (lower_bound, upper_bound)) +``` + +## Optimizing + +Having built a surrogate, we can now use it to search for minima in our original function `f`. + +To optimize using our surrogate we call `surrogate_optimize` method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique. + +```@example linear_surrogate1D +@show surrogate_optimize( + f, SRBF(), lower_bound, upper_bound, my_linear_surr_1D, SobolSample()) +scatter(x, y, label = "Sampled points") +plot!(f, label = "True function", xlims = (lower_bound, upper_bound)) +plot!(my_linear_surr_1D, label = "Surrogate function", xlims = (lower_bound, upper_bound)) +``` + +## Linear Surrogate tutorial (ND) + +First of all we will define the `Egg Holder` function we are going to build a surrogate for. Notice, one how its argument is a vector of numbers, one for each coordinate, and its output is a scalar. + +```@example linear_surrogateND +using Plots # hide +default(c = :matter, legend = false, xlabel = "x", ylabel = "y") # hide +using Surrogates # hide + +function egg(x) + x1 = x[1] + x2 = x[2] + term1 = -(x2 + 47) * sin(sqrt(abs(x2 + x1 / 2 + 47))) + term2 = -x1 * sin(sqrt(abs(x1 - (x2 + 47)))) + 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 `-10, 5`, and `0, 15` for the second dimension. We are taking 50 samples of the space using Sobol Sequences. We then evaluate our function on all of the sampling points. + +```@example linear_surrogateND +n_samples = 50 +lower_bound = [-10.0, 0.0] +upper_bound = [5.0, 15.0] + +xys = sample(n_samples, lower_bound, upper_bound, SobolSample()) +zs = egg.(xys); +``` + +```@example linear_surrogateND +x, y = -10:5, 0:15 # hide +p1 = surface(x, y, (x1, x2) -> egg((x1, x2))) # hide +xs = [xy[1] for xy in xys] # hide +ys = [xy[2] for xy in xys] # hide +scatter!(xs, ys, zs) # hide +p2 = contour(x, y, (x1, x2) -> egg((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 linear_surrogateND +my_linear_ND = LinearSurrogate(xys, zs, lower_bound, upper_bound) +``` + +```@example linear_surrogateND +p1 = surface(x, y, (x, y) -> my_linear_ND([x y])) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide +p2 = contour(x, y, (x, y) -> my_linear_ND([x y])) # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2, title = "Surrogate") # hide +``` + +### Optimizing + +With our surrogate, we can now search for the minima of the function. + +Notice how the new sampled points, which were created during the optimization process, are appended to the `xys` array. +This is why its size changes. + +```@example linear_surrogateND +size(xys) +``` + +```@example linear_surrogateND +surrogate_optimize( + egg, SRBF(), lower_bound, upper_bound, my_linear_ND, SobolSample(), maxiters = 10) +``` + +```@example linear_surrogateND +size(xys) +``` + +```@example linear_surrogateND +p1 = surface(x, y, (x, y) -> my_linear_ND([x y])) # hide +xs = [xy[1] for xy in xys] # hide +ys = [xy[2] for xy in xys] # hide +zs = egg.(xys) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide +p2 = contour(x, y, (x, y) -> my_linear_ND([x y])) # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2) # hide +``` diff --git a/docs/src/Salustowicz.md b/docs/src/Salustowicz.md index 9a4ea0835..130523946 100644 --- a/docs/src/Salustowicz.md +++ b/docs/src/Salustowicz.md @@ -1,64 +1,70 @@ -## 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^2(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 the Salustowicz function with different surrogates: - -```@example salustowicz1D -InverseDistance = InverseDistanceSurrogate(x, y, lower_bound, upper_bound) -lobachevsky_surrogate = LobachevskySurrogate(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, lobachevsky_surrogate.(xs), label="Lobachevsky", 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) -``` +## 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^2(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 the Salustowicz function with different surrogates: + +```@example salustowicz1D +InverseDistance = InverseDistanceSurrogate(x, y, lower_bound, upper_bound) +lobachevsky_surrogate = LobachevskySurrogate( + 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, lobachevsky_surrogate.(xs), label = "Lobachevsky", 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) +``` diff --git a/docs/src/abstractgps.md b/docs/src/abstractgps.md index 2f711af59..f66f27c77 100644 --- a/docs/src/abstractgps.md +++ b/docs/src/abstractgps.md @@ -1,16 +1,20 @@ # Gaussian Process Surrogate Tutorial !!! note - This surrogate requires the 'SurrogatesAbstractGPs' module, which can be added by inputting "]add SurrogatesAbstractGPs" from the Julia command line. + + This surrogate requires the 'SurrogatesAbstractGPs' module, which can be added by inputting "]add SurrogatesAbstractGPs" from the Julia command line. 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. !!! tip + The examples below demonstrate the use of AbstractGPs with out-of-the-box settings without hyperparameter optimization (i.e. without changing parameters like lengthscale, signal variance, and noise variance). Beyond hyperparameter optimization, careful initialization of hyperparameters and priors on the parameters is required for this surrogate to work properly. For more details on how to fit GPs in practice, check out [A Practical Guide to Gaussian Processes](https://infallible-thompson-49de36.netlify.app/). Also see this [example](https://juliagaussianprocesses.github.io/AbstractGPs.jl/stable/examples/1-mauna-loa/#Hyperparameter-Optimization) to understand hyperparameter optimization with AbstractGPs. -## 1D Example -In the example below, the 'gp_surrogate' assignment code can be commented / uncommented to see how the different kernels influence the predictions. + +## 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 @@ -28,13 +32,16 @@ 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) +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 @@ -43,35 +50,37 @@ using Plots using AbstractGPs using SurrogatesAbstractGPs -f(x) = (x-2)^2 +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) +gp_surrogate = AbstractGPSurrogate(x, y) @show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, gp_surrogate, SobolSample()) ``` -Plotting the function and the sampled points: + +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) +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") +default(c = :matter, legend = false, xlabel = "x", ylabel = "y") using Surrogates using AbstractGPs using SurrogatesAbstractGPs - -hypot_func = z -> 3*hypot(z...)+1 +hypot_func = z -> 3 * hypot(z...) + 1 n_samples = 50 lower_bound = [-1.0, -1.0] upper_bound = [1.0, 1.0] @@ -79,35 +88,37 @@ 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))) +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") +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) +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") +scatter!(xs, ys, marker_z = zs) +plot(p1, p2, title = "Surrogate") ``` ```@example abstractgps_tutorialnd -@show gp_surrogate((0.2,0.2)) +@show gp_surrogate((0.2, 0.2)) ``` ```@example abstractgps_tutorialnd -@show hypot_func((0.2,0.2)) +@show hypot_func((0.2, 0.2)) ``` And this is our log marginal posterior predictive probability: + ```@example abstractgps_tutorialnd @show logpdf_surrogate(gp_surrogate) -``` \ No newline at end of file +``` diff --git a/docs/src/ackley.md b/docs/src/ackley.md index fe38d4e9d..e23a6c889 100644 --- a/docs/src/ackley.md +++ b/docs/src/ackley.md @@ -2,7 +2,7 @@ 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`` +Usually the recommended values are: ``a = 20``, ``b = 0.2`` and ``c = 2\pi`` Let's see the 1D case. @@ -16,20 +16,19 @@ Now, let's define the `Ackley` function: ```@example ackley function ackley(x) - a, b, c = 20.0, 0.2, 2.0*π + 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_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) + 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 @@ -37,8 +36,8 @@ 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) +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 @@ -47,21 +46,20 @@ my_loba = LobachevskySurrogate(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_loba.(xs), label="Lobachevsky", legend=:top) - +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_loba.(xs), label = "Lobachevsky", 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,RandomSample()) -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) +surrogate_optimize(ackley, DYCORS(), lb, ub, my_rad, RandomSample()) +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 method successfully finds the minimum. diff --git a/docs/src/cantilever.md b/docs/src/cantilever.md index d259f89dc..c86942971 100644 --- a/docs/src/cantilever.md +++ b/docs/src/cantilever.md @@ -1,9 +1,11 @@ # Cantilever Beam Function + The Cantilever Beam function is defined as: ``f(w,t) = \frac{4L^3}{Ewt}*\sqrt{ (\frac{Y}{t^2})^2 + (\frac{X}{w^2})^2 }`` With parameters L,E,X and Y given. Let's import Surrogates and Plots: + ```@example beam using Surrogates using SurrogatesPolyChaos @@ -12,6 +14,7 @@ default() ``` Define the objective function: + ```@example beam function f(x) t = x[1] @@ -20,56 +23,58 @@ function f(x) E = 2.770674127819261e7 X = 530.8038576066307 Y = 997.8714938733949 - return (4*L^3)/(E*w*t)*sqrt( (Y/t^2)^2 + (X/w^2)^2) + return (4 * L^3) / (E * w * t) * sqrt((Y / t^2)^2 + (X / w^2)^2) end ``` Let's plot it: + ```@example beam n = 100 -lb = [1.0,1.0] -ub = [8.0,8.0] -xys = sample(n,lb,ub,SobolSample()); +lb = [1.0, 1.0] +ub = [8.0, 8.0] +xys = sample(n, lb, ub, SobolSample()); zs = f.(xys); x, y = 0.0:8.0, 0.0:8.0 -p1 = surface(x, y, (x1,x2) -> f((x1,x2))) +p1 = surface(x, y, (x1, x2) -> f((x1, x2))) xs = [xy[1] for xy in xys] ys = [xy[2] for xy in xys] scatter!(xs, ys, zs) # hide -p2 = contour(x, y, (x1,x2) -> f((x1,x2))) +p2 = contour(x, y, (x1, x2) -> f((x1, x2))) scatter!(xs, ys) -plot(p1, p2, title="True function") +plot(p1, p2, title = "True function") ``` - Fitting different surrogates: + ```@example beam -mypoly = PolynomialChaosSurrogate(xys, zs, lb, ub) -loba = LobachevskySurrogate(xys, zs, lb, ub) -rad = RadialBasis(xys,zs,lb,ub) +mypoly = PolynomialChaosSurrogate(xys, zs, lb, ub) +loba = LobachevskySurrogate(xys, zs, lb, ub) +rad = RadialBasis(xys, zs, lb, ub) ``` Plotting: + ```@example beam p1 = surface(x, y, (x, y) -> mypoly([x y])) -scatter!(xs, ys, zs, marker_z=zs) +scatter!(xs, ys, zs, marker_z = zs) p2 = contour(x, y, (x, y) -> mypoly([x y])) -scatter!(xs, ys, marker_z=zs) -plot(p1, p2, title="Polynomial expansion") +scatter!(xs, ys, marker_z = zs) +plot(p1, p2, title = "Polynomial expansion") ``` ```@example beam p1 = surface(x, y, (x, y) -> loba([x y])) -scatter!(xs, ys, zs, marker_z=zs) +scatter!(xs, ys, zs, marker_z = zs) p2 = contour(x, y, (x, y) -> loba([x y])) -scatter!(xs, ys, marker_z=zs) -plot(p1, p2, title="Lobachevsky") +scatter!(xs, ys, marker_z = zs) +plot(p1, p2, title = "Lobachevsky") ``` ```@example beam p1 = surface(x, y, (x, y) -> rad([x y])) -scatter!(xs, ys, zs, marker_z=zs) +scatter!(xs, ys, zs, marker_z = zs) p2 = contour(x, y, (x, y) -> rad([x y])) -scatter!(xs, ys, marker_z=zs) -plot(p1, p2, title="Inverse distance") +scatter!(xs, ys, marker_z = zs) +plot(p1, p2, title = "Inverse distance") ``` diff --git a/docs/src/gek.md b/docs/src/gek.md index 7d780f41d..df8540186 100644 --- a/docs/src/gek.md +++ b/docs/src/gek.md @@ -1,125 +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, which increases proportionally with both the number of inputs and the number of sampling points. - -Let's have a look at 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 and 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); -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 the sampling points. - -```@example GEK_ND -n_samples = 45 -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 -``` +## 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, which increases proportionally with both the number of inputs and the number of sampling points. + +Let's have a look at 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 and 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); +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 the sampling points. + +```@example GEK_ND +n_samples = 45 +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 +``` diff --git a/docs/src/gekpls.md b/docs/src/gekpls.md index 5ea4af75c..f7a474a95 100644 --- a/docs/src/gekpls.md +++ b/docs/src/gekpls.md @@ -1,20 +1,21 @@ ## GEKPLS Surrogate Tutorial -Gradient Enhanced Kriging with Partial Least Squares Method (GEKPLS) is a surrogate modeling technique that brings down computation time and returns improved accuracy for high-dimensional problems. The Julia implementation of GEKPLS is adapted from the Python version by [SMT](https://github.com/SMTorg) which is based on this [paper](https://arxiv.org/pdf/1708.02663.pdf). +Gradient Enhanced Kriging with Partial Least Squares Method (GEKPLS) is a surrogate modeling technique that brings down computation time and returns improved accuracy for high-dimensional problems. The Julia implementation of GEKPLS is adapted from the Python version by [SMT](https://github.com/SMTorg) which is based on this [paper](https://arxiv.org/pdf/1708.02663.pdf). -The following are the inputs when building a GEKPLS surrogate: +The following are the inputs when building a GEKPLS surrogate: -1. x - The vector containing the training points -2. y - The vector containing the training outputs associated with each of the training points -3. grads - The gradients at each of the input X training points -4. n_comp - Number of components to retain for the partial least squares regression (PLS) -5. delta_x - The step size to use for the first order Taylor approximation -6. lb - The lower bound for the training points -7. ub - The upper bound for the training points -8. extra_points - The number of additional points to use for the PLS -9. theta - The hyperparameter to use for the correlation model + 1. x - The vector containing the training points + 2. y - The vector containing the training outputs associated with each of the training points + 3. grads - The gradients at each of the input X training points + 4. n_comp - Number of components to retain for the partial least squares regression (PLS) + 5. delta_x - The step size to use for the first order Taylor approximation + 6. lb - The lower bound for the training points + 7. ub - The upper bound for the training points + 8. extra_points - The number of additional points to use for the PLS + 9. theta - The hyperparameter to use for the correlation model ### Basic GEKPLS Usage + The following example illustrates how to use GEKPLS: ```@example gekpls_water_flow @@ -31,18 +32,19 @@ function water_flow(x) H_l = x[6] L = x[7] K_w = x[8] - log_val = log(r/r_w) - return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l)) + log_val = log(r / r_w) + return (2 * pi * T_u * (H_u - H_l)) / + (log_val * (1 + (2 * L * T_u / (log_val * r_w^2 * K_w)) + T_u / T_l)) end n = 1000 -lb = [0.05,100,63070,990,63.1,700,1120,9855] -ub = [0.15,50000,115600,1110,116,820,1680,12045] -x = sample(n,lb,ub,SobolSample()) +lb = [0.05, 100, 63070, 990, 63.1, 700, 1120, 9855] +ub = [0.15, 50000, 115600, 1110, 116, 820, 1680, 12045] +x = sample(n, lb, ub, SobolSample()) grads = gradient.(water_flow, x) y = water_flow.(x) -n_test = 100 -x_test = sample(n_test,lb,ub,GoldenSample()) +n_test = 100 +x_test = sample(n_test, lb, ub, GoldenSample()) y_true = water_flow.(x_test) n_comp = 2 delta_x = 0.0001 @@ -50,38 +52,37 @@ extra_points = 2 initial_theta = [0.01 for i in 1:n_comp] g = GEKPLS(x, y, grads, n_comp, delta_x, lb, ub, extra_points, initial_theta) y_pred = g.(x_test) -rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) #root mean squared error +rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test)) #root mean squared error println(rmse) #0.0347 - ``` ### Using GEKPLS With Surrogate Optimization + GEKPLS can also be used to find the minimum of a function with the surrogates.jl optimization function. This next example demonstrates how this can be accomplished. ```@example gekpls_optimization - using Surrogates - using Zygote - - function sphere_function(x) - return sum(x .^ 2) - end +using Surrogates +using Zygote - lb = [-5.0, -5.0, -5.0] - ub = [5.0, 5.0, 5.0] - n_comp = 2 - delta_x = 0.0001 - extra_points = 2 - initial_theta = [0.01 for i in 1:n_comp] - n = 100 - x = sample(n, lb, ub, SobolSample()) - grads = gradient.(sphere_function, x) - y = sphere_function.(x) - g = GEKPLS(x, y, grads, n_comp, delta_x, lb, ub, extra_points, initial_theta) - x_point, minima = surrogate_optimize(sphere_function, SRBF(), lb, ub, g, - RandomSample(); maxiters = 20, - num_new_samples = 20, needs_gradient = true) - println(minima) +function sphere_function(x) + return sum(x .^ 2) +end +lb = [-5.0, -5.0, -5.0] +ub = [5.0, 5.0, 5.0] +n_comp = 2 +delta_x = 0.0001 +extra_points = 2 +initial_theta = [0.01 for i in 1:n_comp] +n = 100 +x = sample(n, lb, ub, SobolSample()) +grads = gradient.(sphere_function, x) +y = sphere_function.(x) +g = GEKPLS(x, y, grads, n_comp, delta_x, lb, ub, extra_points, initial_theta) +x_point, minima = surrogate_optimize(sphere_function, SRBF(), lb, ub, g, + RandomSample(); maxiters = 20, + num_new_samples = 20, needs_gradient = true) +println(minima) ``` diff --git a/docs/src/gramacylee.md b/docs/src/gramacylee.md index ad92e00e1..ade423d04 100644 --- a/docs/src/gramacylee.md +++ b/docs/src/gramacylee.md @@ -19,9 +19,9 @@ 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; + term1 = sin(10 * pi * x) / 2 * x + term2 = (x - 1)^4 + y = term1 + term2 end ``` @@ -34,8 +34,9 @@ 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) +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: @@ -44,9 +45,10 @@ Now, let's fit Gramacy & Lee function with different surrogates: my_pol = PolynomialChaosSurrogate(x, y, lower_bound, upper_bound) loba_1 = LobachevskySurrogate(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="Lobachevsky", legend=:top) -plot!(xs, krig.(xs), label="Kriging", legend=:top) +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 = "Lobachevsky", legend = :top) +plot!(xs, krig.(xs), label = "Kriging", legend = :top) ``` diff --git a/docs/src/index.md b/docs/src/index.md index e0484b70d..cc7cb0664 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,35 +1,36 @@ # Surrogates.jl: Surrogate models and optimization for scientific machine learning + A surrogate model is an approximation method that mimics the behavior of a computationally expensive simulation. In more mathematical terms: suppose we are attempting to optimize a function ``\; f(p)``, but each calculation of ``\; f`` is very expensive. It may be the case that we need to solve a PDE for each point or use advanced numerical linear algebra machinery, which is usually costly. The idea is then to develop a surrogate model ``\; g`` which approximates ``\; f`` by training on previous data collected from evaluations of ``\; f``. The construction of a surrogate model can be seen as a three-step process: -1. Sample selection -2. Construction of the surrogate model -3. Surrogate optimization + 1. Sample selection + 2. Construction of the surrogate model + 3. Surrogate optimization The sampling methods are super important for the behavior of the surrogate. Sampling can be done through [QuasiMonteCarlo.jl](https://github.com/SciML/QuasiMonteCarlo.jl), all the functions available there can be used in Surrogates.jl. The available surrogates are: -- Linear -- Radial Basis -- Kriging -- Custom Kriging provided with Stheno -- Neural Network -- Support Vector Machine -- Random Forest -- Second Order Polynomial -- Inverse Distance + - Linear + - Radial Basis + - Kriging + - Custom Kriging provided with Stheno + - Neural Network + - Support Vector Machine + - Random Forest + - Second Order Polynomial + - Inverse Distance After the surrogate is built, we need to optimize it with respect to some objective function. That is, simultaneously looking for a minimum **and** sampling the most unknown region. The available optimization methods are: -- Stochastic RBF (SRBF) -- Lower confidence-bound strategy (LCBS) -- Expected improvement (EI) -- Dynamic coordinate search (DYCORS) + - Stochastic RBF (SRBF) + - Lower confidence-bound strategy (LCBS) + - Expected improvement (EI) + - Dynamic coordinate search (DYCORS) ## Multi-output Surrogates @@ -54,11 +55,11 @@ y = f.(x) Currently, the following are implemented as multi-output surrogates: -- Radial Basis -- Neural Network (via Flux) -- Second Order Polynomial -- Inverse Distance -- Custom Kriging (via Stheno) + - Radial Basis + - Neural Network (via Flux) + - Second Order Polynomial + - Inverse Distance + - Custom Kriging (via Stheno) ## Gradients @@ -67,26 +68,32 @@ of this property, surrogates are useful models for processes which aren't explic differentiable, and can be used as layers in, for instance, Flux models. ## Installation + Surrogates is registered in the Julia General Registry. In the REPL: + ``` using Pkg Pkg.add("Surrogates") ``` + ## Contributing -- Please refer to the - [SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://github.com/SciML/ColPrac/blob/master/README.md) - for guidance on PRs, issues, and other matters relating to contributing to SciML. -- See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions. -- There are a few community forums: - - The #diffeq-bridged and #sciml-bridged channels in the - [Julia Slack](https://julialang.org/slack/) - - The #diffeq-bridged and #sciml-bridged channels in the - [Julia Zulip](https://julialang.zulipchat.com/#narrow/stream/279055-sciml-bridged) - - On the [Julia Discourse forums](https://discourse.julialang.org) - - See also [SciML Community page](https://sciml.ai/community/) + - Please refer to the + [SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://github.com/SciML/ColPrac/blob/master/README.md) + for guidance on PRs, issues, and other matters relating to contributing to SciML. + + - See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions. + - There are a few community forums: + + + The #diffeq-bridged and #sciml-bridged channels in the + [Julia Slack](https://julialang.org/slack/) + + The #diffeq-bridged and #sciml-bridged channels in the + [Julia Zulip](https://julialang.zulipchat.com/#narrow/stream/279055-sciml-bridged) + + On the [Julia Discourse forums](https://discourse.julialang.org) + + See also [SciML Community page](https://sciml.ai/community/) ## Quick example + ```@example using Surrogates num_samples = 10 @@ -94,24 +101,25 @@ lb = 0.0 ub = 10.0 #Sampling -x = sample(num_samples,lb,ub,SobolSample()) -f = x-> log(x)*x^2+x^3 +x = sample(num_samples, lb, ub, SobolSample()) +f = x -> log(x) * x^2 + x^3 y = f.(x) #Creating surrogate alpha = 2.0 n = 6 -my_lobachevsky = LobachevskySurrogate(x,y,lb,ub,alpha=alpha,n=n) +my_lobachevsky = LobachevskySurrogate(x, y, lb, ub, alpha = alpha, n = n) #Approximating value at 5.0 value = my_lobachevsky(5.0) #Adding more data points -surrogate_optimize(f,SRBF(),lb,ub,my_lobachevsky,RandomSample()) +surrogate_optimize(f, SRBF(), lb, ub, my_lobachevsky, RandomSample()) #New approximation value = my_lobachevsky(5.0) ``` + ## Reproducibility ```@raw html diff --git a/docs/src/kriging.md b/docs/src/kriging.md index b6e93fa9b..bc5983c90 100644 --- a/docs/src/kriging.md +++ b/docs/src/kriging.md @@ -5,11 +5,13 @@ Kriging or Gaussian process regression, is a method of interpolation in which th We are going to use a Kriging surrogate to optimize $f(x)=(6x-2)^2sin(12x-4)$. (function from Forrester et al. (2008)). First of all, import `Surrogates` and `Plots`. + ```@example kriging_tutorial1d using Surrogates using Plots default() ``` + ### Sampling We choose to sample f in 4 points between 0 and 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. @@ -28,9 +30,11 @@ xs = lower_bound:0.001:upper_bound x = sample(n_samples, lower_bound, upper_bound, SobolSample()) y = f.(x) -scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound), ylims=(-7, 17)) -plot!(xs, f.(xs), label="True function", legend=:top) +scatter( + x, y, label = "Sampled points", xlims = (lower_bound, upper_bound), ylims = (-7, 17)) +plot!(xs, f.(xs), label = "True function", legend = :top) ``` + ### Building a surrogate With our sampled points, we can build the Kriging surrogate using the `Kriging` function. @@ -40,47 +44,53 @@ With our sampled points, we can build the Kriging surrogate using the `Kriging` ```@example kriging_tutorial1d kriging_surrogate = Kriging(x, y, lower_bound, upper_bound); -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!(xs, kriging_surrogate.(xs), label="Surrogate function", ribbon=p->std_error_at_point(kriging_surrogate, p), legend=:top) +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!(xs, kriging_surrogate.(xs), label = "Surrogate function", + ribbon = p -> std_error_at_point(kriging_surrogate, p), legend = :top) ``` + ### Optimizing + Having built a surrogate, we can now use it to search for minima in our original function `f`. To optimize using our surrogate, we call `surrogate_optimize` method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique. ```@example kriging_tutorial1d -@show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, kriging_surrogate, SobolSample()) +@show surrogate_optimize( + f, SRBF(), lower_bound, upper_bound, kriging_surrogate, SobolSample()) -scatter(x, y, label="Sampled points", ylims=(-7, 7), legend=:top) -plot!(xs, f.(xs), label="True function", legend=:top) -plot!(xs, kriging_surrogate.(xs), label="Surrogate function", ribbon=p->std_error_at_point(kriging_surrogate, p), legend=:top) +scatter(x, y, label = "Sampled points", ylims = (-7, 7), legend = :top) +plot!(xs, f.(xs), label = "True function", legend = :top) +plot!(xs, kriging_surrogate.(xs), label = "Surrogate function", + ribbon = p -> std_error_at_point(kriging_surrogate, p), legend = :top) ``` - ## Kriging surrogate tutorial (ND) First of all, let's define the function we are going to build a surrogate for. Notice how its argument is a vector of numbers, one for each coordinate, and its output is a scalar. ```@example kriging_tutorialnd using Plots # hide -default(c=:matter, legend=false, xlabel="x", ylabel="y") # hide +default(c = :matter, legend = false, xlabel = "x", ylabel = "y") # hide using Surrogates # hide function branin(x) - x1=x[1] - x2=x[2] - a=1; - b=5.1/(4*π^2); - c=5/π; - r=6; - s=10; - t=1/(8π); - a*(x2-b*x1+c*x1-r)^2+s*(1-t)*cos(x1)+s + x1 = x[1] + x2 = x[2] + a = 1 + b = 5.1 / (4 * π^2) + c = 5 / π + r = 6 + s = 10 + t = 1 / (8π) + a * (x2 - b * x1 + c * x1 - r)^2 + s * (1 - t) * cos(x1) + s 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 `-5, 10`, and `0, 15` for the second dimension. We are taking 50 samples of the space using Sobol sequences. We then evaluate our function on all the sampling points. ```@example kriging_tutorialnd @@ -94,30 +104,34 @@ zs = branin.(xys); ```@example kriging_tutorialnd x, y = -5:10, 0:15 # hide -p1 = surface(x, y, (x1,x2) -> branin((x1,x2))) # hide +p1 = surface(x, y, (x1, x2) -> branin((x1, x2))) # hide xs = [xy[1] for xy in xys] # hide ys = [xy[2] for xy in xys] # hide scatter!(xs, ys, zs) # hide -p2 = contour(x, y, (x1,x2) -> branin((x1,x2))) # hide +p2 = contour(x, y, (x1, x2) -> branin((x1, x2))) # hide scatter!(xs, ys) # hide -plot(p1, p2, title="True function") # 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 kriging_tutorialnd -kriging_surrogate = Kriging(xys, zs, lower_bound, upper_bound, p=[2.0, 2.0], theta=[0.03, 0.003]) +kriging_surrogate = Kriging( + xys, zs, lower_bound, upper_bound, p = [2.0, 2.0], theta = [0.03, 0.003]) ``` ```@example kriging_tutorialnd p1 = surface(x, y, (x, y) -> kriging_surrogate([x y])) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide p2 = contour(x, y, (x, y) -> kriging_surrogate([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2, title="Surrogate") # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2, title = "Surrogate") # hide ``` ### Optimizing + With our surrogate, we can now search for the minima of the branin function. Notice how the new sampled points, which were created during the optimization process, are appended to the `xys` array. @@ -126,10 +140,12 @@ This is why its size changes. ```@example kriging_tutorialnd size(xys) ``` -```@example kriging_tutorialnd -surrogate_optimize(branin, SRBF(), lower_bound, upper_bound, kriging_surrogate, SobolSample(); maxiters = 100, num_new_samples = 10) +```@example kriging_tutorialnd +surrogate_optimize(branin, SRBF(), lower_bound, upper_bound, kriging_surrogate, + SobolSample(); maxiters = 100, num_new_samples = 10) ``` + ```@example kriging_tutorialnd size(xys) ``` @@ -139,10 +155,8 @@ p1 = surface(x, y, (x, y) -> kriging_surrogate([x y])) # hide xs = [xy[1] for xy in xys] # hide ys = [xy[2] for xy in xys] # hide zs = branin.(xys) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide p2 = contour(x, y, (x, y) -> kriging_surrogate([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide +scatter!(xs, ys, marker_z = zs) # hide plot(p1, p2) # hide ``` - - diff --git a/docs/src/lobachevsky.md b/docs/src/lobachevsky.md index 76584aa6d..db18b06f7 100644 --- a/docs/src/lobachevsky.md +++ b/docs/src/lobachevsky.md @@ -5,25 +5,28 @@ Lobachevsky splines function is a function that is used for univariate and multi We are going to use a Lobachevsky surrogate to optimize $f(x)=sin(x)+sin(10/3 * x)$. First of all import `Surrogates` and `Plots`. + ```@example LobachevskySurrogate_tutorial using Surrogates using Plots default() ``` + ## Sampling We choose to sample f in 4 points between 0 and 4 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 LobachevskySurrogate_tutorial -f(x) = sin(x) + sin(10/3 * x) +f(x) = sin(x) + sin(10 / 3 * x) n_samples = 5 lower_bound = 1.0 upper_bound = 4.0 x = sample(n_samples, lower_bound, upper_bound, SobolSample()) y = f.(x) -scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound)) -plot!(f, label="True function", xlims=(lower_bound, upper_bound)) +scatter(x, y, label = "Sampled points", xlims = (lower_bound, upper_bound)) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound)) ``` + ## Building a surrogate With our sampled points, we can build the Lobachevsky surrogate using the `LobachevskySurrogate` function. @@ -33,24 +36,30 @@ With our sampled points, we can build the Lobachevsky surrogate using the `Lobac ```@example LobachevskySurrogate_tutorial alpha = 2.0 n = 6 -lobachevsky_surrogate = LobachevskySurrogate(x, y, lower_bound, upper_bound, alpha = 2.0, n = 6) -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound)) -plot!(f, label="True function", xlims=(lower_bound, upper_bound)) -plot!(lobachevsky_surrogate, label="Surrogate function", xlims=(lower_bound, upper_bound)) +lobachevsky_surrogate = LobachevskySurrogate( + x, y, lower_bound, upper_bound, alpha = 2.0, n = 6) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lower_bound, upper_bound)) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound)) +plot!( + lobachevsky_surrogate, label = "Surrogate function", xlims = (lower_bound, upper_bound)) ``` + ## Optimizing + Having built a surrogate, we can now use it to search for minima in our original function `f`. To optimize using our surrogate we call `surrogate_optimize` method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique. ```@example LobachevskySurrogate_tutorial -@show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, lobachevsky_surrogate, SobolSample()) -scatter(x, y, label="Sampled points") -plot!(f, label="True function", xlims=(lower_bound, upper_bound)) -plot!(lobachevsky_surrogate, label="Surrogate function", xlims=(lower_bound, upper_bound)) +@show surrogate_optimize( + f, SRBF(), lower_bound, upper_bound, lobachevsky_surrogate, SobolSample()) +scatter(x, y, label = "Sampled points") +plot!(f, label = "True function", xlims = (lower_bound, upper_bound)) +plot!( + lobachevsky_surrogate, label = "Surrogate function", xlims = (lower_bound, upper_bound)) ``` - In the example below, it shows how to use `lobachevsky_surrogate` for higher dimension problems. # Lobachevsky Surrogate Tutorial (ND): @@ -59,19 +68,18 @@ First of all, we will define the `Schaffer` function we are going to build surro ```@example LobachevskySurrogate_ND using Plots # hide -default(c=:matter, legend=false, xlabel="x", ylabel="y") # hide +default(c = :matter, legend = false, xlabel = "x", ylabel = "y") # hide using Surrogates # hide function schaffer(x) - x1=x[1] - x2=x[2] - fact1 = x1 ^2; - fact2 = x2 ^2; - y = fact1 + fact2; + x1 = x[1] + x2 = x[2] + fact1 = x1^2 + fact2 = x2^2 + y = fact1 + fact2 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, 8`, and `0, 8` for the second dimension. We are taking 60 samples of the space using Sobol Sequences. We then evaluate our function on all of the sampling points. @@ -87,33 +95,34 @@ zs = schaffer.(xys); ```@example LobachevskySurrogate_ND x, y = 0:8, 0:8 # hide -p1 = surface(x, y, (x1,x2) -> schaffer((x1,x2))) # hide +p1 = surface(x, y, (x1, x2) -> schaffer((x1, x2))) # hide xs = [xy[1] for xy in xys] # hide ys = [xy[2] for xy in xys] # hide scatter!(xs, ys, zs) # hide -p2 = contour(x, y, (x1,x2) -> schaffer((x1,x2))) # hide +p2 = contour(x, y, (x1, x2) -> schaffer((x1, x2))) # hide scatter!(xs, ys) # hide -plot(p1, p2, title="True function") # 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 LobachevskySurrogate_ND -Lobachevsky = LobachevskySurrogate(xys, zs, lower_bound, upper_bound, alpha = [2.4,2.4], n=8) +Lobachevsky = LobachevskySurrogate( + xys, zs, lower_bound, upper_bound, alpha = [2.4, 2.4], n = 8) ``` ```@example LobachevskySurrogate_ND p1 = surface(x, y, (x, y) -> Lobachevsky([x y])) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide p2 = contour(x, y, (x, y) -> Lobachevsky([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2, title="Surrogate") # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2, title = "Surrogate") # hide ``` - ## Optimizing + With our surrogate, we can now search for the minima of the function. Notice how the new sampled points, which were created during the optimization process, are appended to the `xys` array. @@ -122,9 +131,12 @@ This is why its size changes. ```@example LobachevskySurrogate_ND size(Lobachevsky.x) ``` + ```@example LobachevskySurrogate_ND -surrogate_optimize(schaffer, SRBF(), lower_bound, upper_bound, Lobachevsky, SobolSample(), maxiters=1, num_new_samples=10) +surrogate_optimize(schaffer, SRBF(), lower_bound, upper_bound, Lobachevsky, + SobolSample(), maxiters = 1, num_new_samples = 10) ``` + ```@example LobachevskySurrogate_ND size(Lobachevsky.x) ``` @@ -135,8 +147,8 @@ xys = Lobachevsky.x # hide xs = [i[1] for i in xys] # hide ys = [i[2] for i in xys] # hide zs = schaffer.(xys) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide p2 = contour(x, y, (x, y) -> Lobachevsky([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide +scatter!(xs, ys, marker_z = zs) # hide plot(p1, p2) # hide ``` diff --git a/docs/src/lp.md b/docs/src/lp.md index 1e47e97da..6cbd49153 100644 --- a/docs/src/lp.md +++ b/docs/src/lp.md @@ -1,9 +1,10 @@ # Lp norm function + The Lp norm function is defined as: ``f(x) = \sqrt[p]{ \sum_{i=1}^d \vert x_i \vert ^p}`` - Let's import Surrogates and Plots: + ```@example lp using Surrogates using SurrogatesPolyChaos @@ -13,33 +14,38 @@ default() ``` Define the objective function: + ```@example lp -function f(x,p) - return norm(x,p) +function f(x, p) + return norm(x, p) end ``` Let's see a simple 1D case: + ```@example lp n = 30 lb = -5.0 ub = 5.0 p = 1.3 -x = sample(n,lb,ub,SobolSample()) -y = f.(x,p) +x = sample(n, lb, ub, SobolSample()) +y = f.(x, p) xs = lb:0.001:ub -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lb, ub), ylims=(0, 5), legend=:top) -plot!(xs,f.(xs,p), label="True function", legend=:top) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lb, ub), ylims = (0, 5), legend = :top) +plot!(xs, f.(xs, p), label = "True function", legend = :top) ``` Fitting different surrogates: + ```@example lp -my_pol = PolynomialChaosSurrogate(x,y,lb,ub) -loba_1 = LobachevskySurrogate(x,y,lb,ub) -krig = Kriging(x,y,lb,ub) -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lb, ub), ylims=(0, 5), legend=:top) -plot!(xs,f.(xs,p), label="True function", legend=:top) -plot!(xs, my_pol.(xs), label="Polynomial expansion", legend=:top) -plot!(xs, loba_1.(xs), label="Lobachevsky", legend=:top) -plot!(xs, krig.(xs), label="Kriging", legend=:top) +my_pol = PolynomialChaosSurrogate(x, y, lb, ub) +loba_1 = LobachevskySurrogate(x, y, lb, ub) +krig = Kriging(x, y, lb, ub) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lb, ub), ylims = (0, 5), legend = :top) +plot!(xs, f.(xs, p), label = "True function", legend = :top) +plot!(xs, my_pol.(xs), label = "Polynomial expansion", legend = :top) +plot!(xs, loba_1.(xs), label = "Lobachevsky", legend = :top) +plot!(xs, krig.(xs), label = "Kriging", legend = :top) ``` diff --git a/docs/src/moe.md b/docs/src/moe.md index a448bc6e3..42f27fa38 100644 --- a/docs/src/moe.md +++ b/docs/src/moe.md @@ -1,7 +1,8 @@ ## Mixture of Experts (MOE) !!! note - This surrogate requires the 'SurrogatesMOE' module, which can be added by inputting "]add SurrogatesMOE" from the Julia command line. + + This surrogate requires the 'SurrogatesMOE' module, which can be added by inputting "]add SurrogatesMOE" from the Julia command line. The Mixture of Experts (MOE) Surrogate model represents the interpolating function as a combination of other surrogate models. SurrogatesMOE is a Julia implementation of the [Python version from SMT](https://smt.readthedocs.io/en/latest/_src_docs/applications/moe.html). @@ -33,7 +34,8 @@ lb = -1.0 ub = 1.0 x = sample(50, lb, ub, SobolSample()) y = discont_1D.(x) -scatter(x, y, label="Sampled Points", xlims=(lb, ub), ylims=(-6.0, 7.0), legend=:top) +scatter( + x, y, label = "Sampled Points", xlims = (lb, ub), ylims = (-6.0, 7.0), legend = :top) ``` How does a regular surrogate perform on such a dataset? @@ -47,17 +49,17 @@ As we can see, the prediction is far from the ground truth. Now, how does the MO ```@example MOE_1D expert_types = [ - RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, - sparse = false), - RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, - sparse = false) - ] + RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, + sparse = false), + RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, + sparse = false) +] MOE_1D_RAD_RAD = MOE(x, y, expert_types) MOE_at0 = MOE_1D_RAD_RAD(0.0) ``` -As we can see, the accuracy is significantly better. +As we can see, the accuracy is significantly better. ### Under the Hood - How SurrogatesMOE Works @@ -98,9 +100,9 @@ x_test = sample(10, lb, ub, GoldenSample()) expert_types = [ RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, - sparse = false), + sparse = false), RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, - sparse = false), + sparse = false) ] moe_nd_rad_rad = MOE(x, y, expert_types, ndim = 2) moe_pred_vals = moe_nd_rad_rad.(x_test) @@ -127,9 +129,9 @@ expert_types = [ #With 3 Surrogates expert_types = [ RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, - sparse = false), + sparse = false), LinearStructure(), - InverseDistanceStructure(p = 1.0), + InverseDistanceStructure(p = 1.0) ] nothing # hide ``` diff --git a/docs/src/multi_objective_opt.md b/docs/src/multi_objective_opt.md index eb066664d..e4dfb2bf1 100644 --- a/docs/src/multi_objective_opt.md +++ b/docs/src/multi_objective_opt.md @@ -1,4 +1,4 @@ -# Multi-objective optimization +# Multi-objective optimization ## Case 1: Non-colliding objective functions @@ -6,42 +6,44 @@ using Surrogates m = 10 -f = x -> [x^i for i = 1:m] +f = x -> [x^i for i in 1:m] lb = 1.0 ub = 10.0 -x = sample(50, lb, ub, GoldenSample()) -y = f.(x) +x = sample(50, lb, ub, GoldenSample()) +y = f.(x) my_radial_basis_ego = RadialBasis(x, y, lb, ub) -pareto_set, pareto_front = surrogate_optimize(f, SMB(),lb,ub,my_radial_basis_ego,SobolSample(); maxiters = 10, n_new_look = 100) +pareto_set, pareto_front = surrogate_optimize( + f, SMB(), lb, ub, my_radial_basis_ego, SobolSample(); maxiters = 10, n_new_look = 100) m = 5 -f = x -> [x^i for i =1:m] +f = x -> [x^i for i in 1:m] lb = 1.0 ub = 10.0 -x = sample(50, lb, ub, SobolSample()) -y = f.(x) +x = sample(50, lb, ub, SobolSample()) +y = f.(x) my_radial_basis_rtea = RadialBasis(x, y, lb, ub) Z = 0.8 K = 2 p_cross = 0.5 n_c = 1.0 sigma = 1.5 -surrogate_optimize(f,RTEA(Z,K,p_cross,n_c,sigma),lb,ub,my_radial_basis_rtea,SobolSample()) - +surrogate_optimize( + f, RTEA(Z, K, p_cross, n_c, sigma), lb, ub, my_radial_basis_rtea, SobolSample()) ``` ## Case 2: objective functions with conflicting minima ```@example multi_obj -f = x -> [sqrt((x[1] - 4)^2 + 25*(x[2])^2), - sqrt((x[1]+4)^2 + 25*(x[2])^2), - sqrt((x[1]-3)^2 + (x[2]-1)^2)] -lb = [2.5,-0.5] -ub = [3.5,0.5] -x = sample(50, lb, ub, SobolSample()) -y = f.(x) +f = x -> [sqrt((x[1] - 4)^2 + 25 * (x[2])^2), + sqrt((x[1] + 4)^2 + 25 * (x[2])^2), + sqrt((x[1] - 3)^2 + (x[2] - 1)^2)] +lb = [2.5, -0.5] +ub = [3.5, 0.5] +x = sample(50, lb, ub, SobolSample()) +y = f.(x) my_radial_basis_ego = RadialBasis(x, y, lb, ub) #I can find my pareto set and pareto front by calling again the surrogate_optimize function: -pareto_set, pareto_front = surrogate_optimize(f,SMB(),lb,ub,my_radial_basis_ego, SobolSample(); maxiters = 10, n_new_look = 100); +pareto_set, pareto_front = surrogate_optimize( + f, SMB(), lb, ub, my_radial_basis_ego, SobolSample(); maxiters = 10, n_new_look = 100); ``` diff --git a/docs/src/neural.md b/docs/src/neural.md index 3a8311b8c..fe3d2417a 100644 --- a/docs/src/neural.md +++ b/docs/src/neural.md @@ -1,6 +1,8 @@ # Neural network tutorial + !!! note - This surrogate requires the 'SurrogatesFlux' module, which can be added by inputting "]add SurrogatesFlux" from the Julia command line. + + This surrogate requires the 'SurrogatesFlux' module, which can be added by inputting "]add SurrogatesFlux" from the Julia command line. It's possible to define a neural network as a surrogate, using Flux. This is useful because we can call optimization methods on it. @@ -9,21 +11,20 @@ First of all we will define the `Schaffer` function we are going to build a surr ```@example Neural_surrogate using Plots -default(c=:matter, legend=false, xlabel="x", ylabel="y") # hide +default(c = :matter, legend = false, xlabel = "x", ylabel = "y") # hide using Surrogates using Flux using SurrogatesFlux function schaffer(x) - x1=x[1] - x2=x[2] - fact1 = x1 ^2; - fact2 = x2 ^2; - y = fact1 + fact2; + x1 = x[1] + x2 = x[2] + fact1 = x1^2 + fact2 = x2^2 + y = fact1 + fact2 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, 8`, and `0, 8` for the second dimension. We are taking 60 samples of the space using Sobol Sequences. We then evaluate our function on all the sampling points. @@ -39,31 +40,34 @@ zs = schaffer.(xys); ```@example Neural_surrogate x, y = 0:8, 0:8 # hide -p1 = surface(x, y, (x1,x2) -> schaffer((x1,x2))) # hide +p1 = surface(x, y, (x1, x2) -> schaffer((x1, x2))) # hide xs = [xy[1] for xy in xys] # hide ys = [xy[2] for xy in xys] # hide scatter!(xs, ys, zs) # hide -p2 = contour(x, y, (x1,x2) -> schaffer((x1,x2))) # hide +p2 = contour(x, y, (x1, x2) -> schaffer((x1, x2))) # hide scatter!(xs, ys) # hide -plot(p1, p2, title="True function") # hide +plot(p1, p2, title = "True function") # hide ``` - ## Building a surrogate + You can specify your own model, optimization function, loss functions, and epochs. As always, getting the model right is the hardest thing. ```@example Neural_surrogate model1 = Chain( - Dense(2, 5, σ), - Dense(5,2,σ), - Dense(2, 1) + Dense(2, 5, σ), + Dense(5, 2, σ), + Dense(2, 1) ) neural = NeuralSurrogate(xys, zs, lower_bound, upper_bound, model = model1, n_echos = 10) ``` ## Optimization + We can now call an optimization function on the neural network: + ```@example Neural_surrogate -surrogate_optimize(schaffer, SRBF(), lower_bound, upper_bound, neural, SobolSample(), maxiters=20, num_new_samples=10) +surrogate_optimize(schaffer, SRBF(), lower_bound, upper_bound, neural, + SobolSample(), maxiters = 20, num_new_samples = 10) ``` diff --git a/docs/src/optimizations.md b/docs/src/optimizations.md index 90b517867..39acfe4c8 100644 --- a/docs/src/optimizations.md +++ b/docs/src/optimizations.md @@ -1,32 +1,40 @@ # Optimization techniques -* SRBF + - SRBF + ```@docs surrogate_optimize(obj::Function,::SRBF,lb,ub,surr::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100) ``` -* LCBS + - LCBS + ```@docs surrogate_optimize(obj::Function,::LCBS,lb,ub,krig,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100) ``` -* EI + - EI + ```@docs surrogate_optimize(obj::Function,::EI,lb,ub,krig,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100) ``` -* DYCORS + - DYCORS + ```@docs surrogate_optimize(obj::Function,::DYCORS,lb,ub,surrn::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100) ``` -* SOP + - SOP + ```@docs surrogate_optimize(obj::Function,sop1::SOP,lb::Number,ub::Number,surrSOP::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=min(500*1,5000)) ``` + ## Adding another optimization method + To add another optimization method, you just need to define a new SurrogateOptimizationAlgorithm and write its corresponding algorithm, overloading the following: + ``` surrogate_optimize(obj::Function,::NewOptimizationType,lb,ub,surr::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100) ``` diff --git a/docs/src/parallel.md b/docs/src/parallel.md index fccffbbdd..39383e501 100755 --- a/docs/src/parallel.md +++ b/docs/src/parallel.md @@ -1,6 +1,6 @@ # Parallel Optimization -There are some situations where it can be beneficial to run multiple optimizations in parallel. For example, if your objective function is very expensive to evaluate, you may want to run multiple evaluations in parallel. +There are some situations where it can be beneficial to run multiple optimizations in parallel. For example, if your objective function is very expensive to evaluate, you may want to run multiple evaluations in parallel. ## Ask-Tell Interface @@ -9,30 +9,38 @@ To enable parallel optimization, we make use of an Ask-Tell interface. The user ## Virtual Points To ensure that points of interest returned by `potential_optimal_points` are sufficiently far from each other, the function makes use of *virtual points*. They are used as follows: -1. `potential_optimal_points` is told to return `n` points. -2. The point with the highest merit function value is selected. -3. This point is now treated as a virtual point and is assigned a temporary value that changes the landscape of the merit function. How the temporary value is chosen depends on the strategy used. (see below) -4. The point with the new highest merit is selected. -5. The process is repeated until `n` points have been selected. + + 1. `potential_optimal_points` is told to return `n` points. + 2. The point with the highest merit function value is selected. + 3. This point is now treated as a virtual point and is assigned a temporary value that changes the landscape of the merit function. How the temporary value is chosen depends on the strategy used. (see below) + 4. The point with the new highest merit is selected. + 5. The process is repeated until `n` points have been selected. The following strategies are available for virtual point selection for all optimization algorithms: -- "Minimum Constant Liar (MinimumConstantLiar)": - - The virtual point is assigned using the lowest known value of the merit function across all evaluated points. -- "Mean Constant Liar (MeanConstantLiar)": - - The virtual point is assigned using the mean of the merit function across all evaluated points. -- "Maximum Constant Liar (MaximumConstantLiar)": - - The virtual point is assigned using the greatest known value of the merit function across all evaluated points. + - "Minimum Constant Liar (MinimumConstantLiar)": + + + The virtual point is assigned using the lowest known value of the merit function across all evaluated points. + + - "Mean Constant Liar (MeanConstantLiar)": + + + The virtual point is assigned using the mean of the merit function across all evaluated points. + - "Maximum Constant Liar (MaximumConstantLiar)": + + + The virtual point is assigned using the greatest known value of the merit function across all evaluated points. -For Kriging surrogates, specifically, the above and following strategies are available: +For Kriging surrogates, specifically, the above and following strategies are available: -- "Kriging Believer (KrigingBeliever): - - The virtual point is assigned using the mean of the Kriging surrogate at the virtual point. -- "Kriging Believer Upper Bound (KrigingBelieverUpperBound)": - - The virtual point is assigned using 3$\sigma$ above the mean of the Kriging surrogate at the virtual point. -- "Kriging Believer Lower Bound (KrigingBelieverLowerBound)": - - The virtual point is assigned using 3$\sigma$ below the mean of the Kriging surrogate at the virtual point. + - "Kriging Believer (KrigingBeliever): + + + The virtual point is assigned using the mean of the Kriging surrogate at the virtual point. + - "Kriging Believer Upper Bound (KrigingBelieverUpperBound)": + + + The virtual point is assigned using 3$\sigma$ above the mean of the Kriging surrogate at the virtual point. + - "Kriging Believer Lower Bound (KrigingBelieverLowerBound)": + + + The virtual point is assigned using 3$\sigma$ below the mean of the Kriging surrogate at the virtual point. In general, MinimumConstantLiar and KrigingBelieverLowerBound tend to favor exploitation, while MaximumConstantLiar and KrigingBelieverUpperBound tend to favor exploration. MeanConstantLiar and KrigingBeliever tend to be compromises between the two. @@ -50,7 +58,8 @@ y = f.(x) my_k = Kriging(x, y, lb, ub) for _ in 1:10 - new_x, eis = potential_optimal_points(EI(), MeanConstantLiar(), lb, ub, my_k, SobolSample(), 3) + new_x, eis = potential_optimal_points( + EI(), MeanConstantLiar(), lb, ub, my_k, SobolSample(), 3) add_point!(my_k, new_x, f.(new_x)) end ``` diff --git a/docs/src/polychaos.md b/docs/src/polychaos.md index 516445092..10edde266 100644 --- a/docs/src/polychaos.md +++ b/docs/src/polychaos.md @@ -1,12 +1,14 @@ # Polynomial chaos surrogate !!! note - This surrogate requires the 'SurrogatesPolyChaos' module which can be added by inputting "]add SurrogatesPolyChaos" from the Julia command line. + + This surrogate requires the 'SurrogatesPolyChaos' module which can be added by inputting "]add SurrogatesPolyChaos" from the Julia command line. We can create a surrogate using a polynomial expansion, with a different polynomial basis depending on the distribution of the data we are trying to fit. Under the hood, PolyChaos.jl has been used. It is possible to specify a type of polynomial for each dimension of the problem. + ### Sampling We choose to sample f in 25 points between 0 and 10 using the `sample` function. The sampling points are chosen using a Low Discrepancy. This can be done by passing `HaltonSample()` to the `sample` function. @@ -20,21 +22,22 @@ default() n = 20 lower_bound = 1.0 upper_bound = 6.0 -x = sample(n,lower_bound,upper_bound,HaltonSample()) -f = x -> log(x)*x + sin(x) +x = sample(n, lower_bound, upper_bound, HaltonSample()) +f = x -> log(x) * x + sin(x) y = f.(x) -scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) +scatter(x, y, 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 ```@example polychaos -poly1 = PolynomialChaosSurrogate(x,y,lower_bound,upper_bound) -poly2 = PolynomialChaosSurrogate(x,y,lower_bound,upper_bound, op = SurrogatesPolyChaos.GaussOrthoPoly(5)) -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) -plot!(poly1, label="First polynomial", xlims=(lower_bound, upper_bound), legend=:top) -plot!(poly2, label="Second polynomial", xlims=(lower_bound, upper_bound), legend=:top) +poly1 = PolynomialChaosSurrogate(x, y, lower_bound, upper_bound) +poly2 = PolynomialChaosSurrogate( + x, y, lower_bound, upper_bound, op = SurrogatesPolyChaos.GaussOrthoPoly(5)) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lower_bound, upper_bound), legend = :top) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) +plot!(poly1, label = "First polynomial", xlims = (lower_bound, upper_bound), legend = :top) +plot!(poly2, label = "Second polynomial", xlims = (lower_bound, upper_bound), legend = :top) ``` diff --git a/docs/src/radials.md b/docs/src/radials.md index 420f6683c..0854ed6e0 100644 --- a/docs/src/radials.md +++ b/docs/src/radials.md @@ -1,159 +1,166 @@ -## Radial Surrogates -The Radial Basis Surrogate model represents the interpolating function as a linear combination of basis functions, one for each training point. Let's start with something easy to get our hands dirty. I want to build a surrogate for: - -```math -f(x) = \log(x) \cdot x^2+x^3 -``` - -Let's choose the Radial Basis Surrogate for 1D. First of all we have to import these two packages: `Surrogates` and `Plots`, - -```@example RadialBasisSurrogate -using Surrogates -using Plots -default() -``` - -We choose to sample f in 30 points between 5 to 25 using `sample` function. The sampling points are chosen using a Sobol sequence, this can be done by passing `SobolSample()` to the `sample` function. - -```@example RadialBasisSurrogate -f(x) = log(x)*x^2 + x^3 -n_samples = 30 -lower_bound = 5.0 -upper_bound = 25.0 -x = sort(sample(n_samples, lower_bound, upper_bound, SobolSample())) -y = f.(x) -scatter(x, y, label="Sampled Points", xlims=(lower_bound, upper_bound), legend=:top) -plot!(x, y, label="True function", legend=:top) -``` - - -## Building Surrogate - -With our sampled points we can build the **Radial Surrogate** using the `RadialBasis` function. - -We can simply calculate `radial_surrogate` for any value. - -```@example RadialBasisSurrogate -radial_surrogate = RadialBasis(x, y, lower_bound, upper_bound) -val = radial_surrogate(5.4) -``` - -We can also use cubic radial basis functions. - -```@example RadialBasisSurrogate -radial_surrogate = RadialBasis(x, y, lower_bound, upper_bound, rad = cubicRadial()) -val = radial_surrogate(5.4) -``` - -Currently, available radial basis functions are `linearRadial` (the default), `cubicRadial`, `multiquadricRadial`, and `thinplateRadial`. - -Now, we will simply plot `radial_surrogate`: - -```@example RadialBasisSurrogate -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) -plot!(radial_surrogate, label="Surrogate function", xlims=(lower_bound, upper_bound), legend=:top) -``` - - -## Optimizing - -Having built a surrogate, we can now use it to search for minima in our original function `f`. - -To optimize using our surrogate, we call `surrogate_optimize` method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique. - -```@example RadialBasisSurrogate -@show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, radial_surrogate, SobolSample()) -scatter(x, y, label="Sampled points", legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) -plot!(radial_surrogate, label="Surrogate function", xlims=(lower_bound, upper_bound), legend=:top) -``` - - -## Radial Basis Surrogate tutorial (ND) - -First of all, we will define the `Booth` function we are going to build the surrogate for: - -$f(x) = (x_1 + 2*x_2 - 7)^2 + (2*x_1 + x_2 - 5)^2$ - - Notice, how its argument is a vector of numbers, one for each coordinate, and its output is a scalar. - -```@example RadialBasisSurrogateND -using Plots # hide -default(c=:matter, legend=false, xlabel="x", ylabel="y") # hide -using Surrogates # hide - -function booth(x) - x1=x[1] - x2=x[2] - term1 = (x1 + 2*x2 - 7)^2; - term2 = (2*x1 + x2 - 5)^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 `-5, 10`, and `0, 15` 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 RadialBasisSurrogateND -n_samples = 80 -lower_bound = [-5.0, 0.0] -upper_bound = [10.0, 15.0] - -xys = sample(n_samples, lower_bound, upper_bound, SobolSample()) -zs = booth.(xys); -``` - -```@example RadialBasisSurrogateND -x, y = -5.0:10.0, 0.0:15.0 # hide -p1 = surface(x, y, (x1,x2) -> booth((x1,x2))) # hide -xs = [xy[1] for xy in xys] # hide -ys = [xy[2] for xy in xys] # hide -scatter!(xs, ys, zs) # hide -p2 = contour(x, y, (x1,x2) -> booth((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 RadialBasisSurrogateND -radial_basis = RadialBasis(xys, zs, lower_bound, upper_bound) -``` - -```@example RadialBasisSurrogateND -p1 = surface(x, y, (x, y) -> radial_basis([x y])) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide -p2 = contour(x, y, (x, y) -> radial_basis([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2, title="Surrogate") # hide -``` - -### Optimizing -With our surrogate, we can now search for the minima of the function. - -Notice how the new sampled points, which were created during the optimization process, are appended to the `xys` array. -This is why its size changes. - -```@example RadialBasisSurrogateND -size(xys) -``` -```@example RadialBasisSurrogateND -surrogate_optimize(booth, SRBF(), lower_bound, upper_bound, radial_basis, RandomSample(), maxiters=50) -``` -```@example RadialBasisSurrogateND -size(xys) -``` - -```@example RadialBasisSurrogateND -p1 = surface(x, y, (x, y) -> radial_basis([x y])) # hide -xs = [xy[1] for xy in xys] # hide -ys = [xy[2] for xy in xys] # hide -zs = booth.(xys) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide -p2 = contour(x, y, (x, y) -> radial_basis([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2) # hide -``` +## Radial Surrogates + +The Radial Basis Surrogate model represents the interpolating function as a linear combination of basis functions, one for each training point. Let's start with something easy to get our hands dirty. I want to build a surrogate for: + +```math +f(x) = \log(x) \cdot x^2+x^3 +``` + +Let's choose the Radial Basis Surrogate for 1D. First of all we have to import these two packages: `Surrogates` and `Plots`, + +```@example RadialBasisSurrogate +using Surrogates +using Plots +default() +``` + +We choose to sample f in 30 points between 5 to 25 using `sample` function. The sampling points are chosen using a Sobol sequence, this can be done by passing `SobolSample()` to the `sample` function. + +```@example RadialBasisSurrogate +f(x) = log(x) * x^2 + x^3 +n_samples = 30 +lower_bound = 5.0 +upper_bound = 25.0 +x = sort(sample(n_samples, lower_bound, upper_bound, SobolSample())) +y = f.(x) +scatter(x, y, label = "Sampled Points", xlims = (lower_bound, upper_bound), legend = :top) +plot!(x, y, label = "True function", legend = :top) +``` + +## Building Surrogate + +With our sampled points we can build the **Radial Surrogate** using the `RadialBasis` function. + +We can simply calculate `radial_surrogate` for any value. + +```@example RadialBasisSurrogate +radial_surrogate = RadialBasis(x, y, lower_bound, upper_bound) +val = radial_surrogate(5.4) +``` + +We can also use cubic radial basis functions. + +```@example RadialBasisSurrogate +radial_surrogate = RadialBasis(x, y, lower_bound, upper_bound, rad = cubicRadial()) +val = radial_surrogate(5.4) +``` + +Currently, available radial basis functions are `linearRadial` (the default), `cubicRadial`, `multiquadricRadial`, and `thinplateRadial`. + +Now, we will simply plot `radial_surrogate`: + +```@example RadialBasisSurrogate +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lower_bound, upper_bound), legend = :top) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) +plot!(radial_surrogate, label = "Surrogate function", + xlims = (lower_bound, upper_bound), legend = :top) +``` + +## Optimizing + +Having built a surrogate, we can now use it to search for minima in our original function `f`. + +To optimize using our surrogate, we call `surrogate_optimize` method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique. + +```@example RadialBasisSurrogate +@show surrogate_optimize( + f, SRBF(), lower_bound, upper_bound, radial_surrogate, SobolSample()) +scatter(x, y, label = "Sampled points", legend = :top) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) +plot!(radial_surrogate, label = "Surrogate function", + xlims = (lower_bound, upper_bound), legend = :top) +``` + +## Radial Basis Surrogate tutorial (ND) + +First of all, we will define the `Booth` function we are going to build the surrogate for: + +$f(x) = (x_1 + 2*x_2 - 7)^2 + (2*x_1 + x_2 - 5)^2$ + +Notice, how its argument is a vector of numbers, one for each coordinate, and its output is a scalar. + +```@example RadialBasisSurrogateND +using Plots # hide +default(c = :matter, legend = false, xlabel = "x", ylabel = "y") # hide +using Surrogates # hide + +function booth(x) + x1 = x[1] + x2 = x[2] + term1 = (x1 + 2 * x2 - 7)^2 + term2 = (2 * x1 + x2 - 5)^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 `-5, 10`, and `0, 15` 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 RadialBasisSurrogateND +n_samples = 80 +lower_bound = [-5.0, 0.0] +upper_bound = [10.0, 15.0] + +xys = sample(n_samples, lower_bound, upper_bound, SobolSample()) +zs = booth.(xys); +``` + +```@example RadialBasisSurrogateND +x, y = -5.0:10.0, 0.0:15.0 # hide +p1 = surface(x, y, (x1, x2) -> booth((x1, x2))) # hide +xs = [xy[1] for xy in xys] # hide +ys = [xy[2] for xy in xys] # hide +scatter!(xs, ys, zs) # hide +p2 = contour(x, y, (x1, x2) -> booth((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 RadialBasisSurrogateND +radial_basis = RadialBasis(xys, zs, lower_bound, upper_bound) +``` + +```@example RadialBasisSurrogateND +p1 = surface(x, y, (x, y) -> radial_basis([x y])) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide +p2 = contour(x, y, (x, y) -> radial_basis([x y])) # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2, title = "Surrogate") # hide +``` + +### Optimizing + +With our surrogate, we can now search for the minima of the function. + +Notice how the new sampled points, which were created during the optimization process, are appended to the `xys` array. +This is why its size changes. + +```@example RadialBasisSurrogateND +size(xys) +``` + +```@example RadialBasisSurrogateND +surrogate_optimize( + booth, SRBF(), lower_bound, upper_bound, radial_basis, RandomSample(), maxiters = 50) +``` + +```@example RadialBasisSurrogateND +size(xys) +``` + +```@example RadialBasisSurrogateND +p1 = surface(x, y, (x, y) -> radial_basis([x y])) # hide +xs = [xy[1] for xy in xys] # hide +ys = [xy[2] for xy in xys] # hide +zs = booth.(xys) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide +p2 = contour(x, y, (x, y) -> radial_basis([x y])) # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2) # hide +``` diff --git a/docs/src/randomforest.md b/docs/src/randomforest.md index 4086a0b51..e56584745 100644 --- a/docs/src/randomforest.md +++ b/docs/src/randomforest.md @@ -1,33 +1,37 @@ ## Random forests surrogate tutorial !!! note - This surrogate requires the 'SurrogatesRandomForest' module, which can be added by inputting "]add SurrogatesRandomForest" from the Julia command line. + + This surrogate requires the 'SurrogatesRandomForest' module, which can be added by inputting "]add SurrogatesRandomForest" from the Julia command line. Random forests is a supervised learning algorithm that randomly creates and merges multiple decision trees into one forest. We are going to use a random forests surrogate to optimize $f(x)=sin(x)+sin(10/3 * x)$. First of all import `Surrogates` and `Plots`. + ```@example RandomForestSurrogate_tutorial using Surrogates using SurrogatesRandomForest using Plots default() ``` + ### Sampling We choose to sample f in 4 points between 0 and 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 RandomForestSurrogate_tutorial -f(x) = sin(x) + sin(10/3 * x) +f(x) = sin(x) + sin(10 / 3 * x) n_samples = 5 lower_bound = 2.7 upper_bound = 7.5 x = sample(n_samples, lower_bound, upper_bound, SobolSample()) y = f.(x) -scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound)) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) +scatter(x, y, label = "Sampled points", xlims = (lower_bound, upper_bound)) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) ``` + ### Building a surrogate With our sampled points, we can build the Random forests surrogate using the `RandomForestSurrogate` function. @@ -37,43 +41,48 @@ using the parameter num_round ```@example RandomForestSurrogate_tutorial num_round = 2 -randomforest_surrogate = RandomForestSurrogate(x ,y ,lower_bound, upper_bound, num_round = 2) -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) -plot!(randomforest_surrogate, label="Surrogate function", xlims=(lower_bound, upper_bound), legend=:top) +randomforest_surrogate = RandomForestSurrogate( + x, y, lower_bound, upper_bound, num_round = 2) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lower_bound, upper_bound), legend = :top) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) +plot!(randomforest_surrogate, label = "Surrogate function", + xlims = (lower_bound, upper_bound), legend = :top) ``` + ### Optimizing + Having built a surrogate, we can now use it to search for minima in our original function `f`. To optimize using our surrogate, we call `surrogate_optimize` method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique. ```@example RandomForestSurrogate_tutorial -@show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, randomforest_surrogate, SobolSample()) -scatter(x, y, label="Sampled points") -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) -plot!(randomforest_surrogate, label="Surrogate function", xlims=(lower_bound, upper_bound), legend=:top) +@show surrogate_optimize( + f, SRBF(), lower_bound, upper_bound, randomforest_surrogate, SobolSample()) +scatter(x, y, label = "Sampled points") +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) +plot!(randomforest_surrogate, label = "Surrogate function", + xlims = (lower_bound, upper_bound), legend = :top) ``` - ## Random Forest ND First of all we will define the `Bukin Function N. 6` function we are going to build a surrogate for. ```@example RandomForestSurrogateND using Plots # hide -default(c=:matter, legend=false, xlabel="x", ylabel="y") # hide +default(c = :matter, legend = false, xlabel = "x", ylabel = "y") # hide using Surrogates # hide function bukin6(x) - x1=x[1] - x2=x[2] - term1 = 100 * sqrt(abs(x2 - 0.01*x1^2)); - term2 = 0.01 * abs(x1+10); - y = term1 + term2; + x1 = x[1] + x2 = x[2] + term1 = 100 * sqrt(abs(x2 - 0.01 * x1^2)) + term2 = 0.01 * abs(x1 + 10) + 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 `-5, 10`, and `0, 15` for the second dimension. We are taking 50 samples of the space using Sobol Sequences. We then evaluate our function on all the sampling points. @@ -89,34 +98,34 @@ zs = bukin6.(xys); ```@example RandomForestSurrogateND x, y = -5:10, 0:15 # hide -p1 = surface(x, y, (x1,x2) -> bukin6((x1,x2))) # hide +p1 = surface(x, y, (x1, x2) -> bukin6((x1, x2))) # hide xs = [xy[1] for xy in xys] # hide ys = [xy[2] for xy in xys] # hide scatter!(xs, ys, zs) # hide -p2 = contour(x, y, (x1,x2) -> bukin6((x1,x2))) # hide +p2 = contour(x, y, (x1, x2) -> bukin6((x1, x2))) # hide scatter!(xs, ys) # hide -plot(p1, p2, title="True function") # 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 RandomForestSurrogateND using SurrogatesRandomForest -RandomForest = RandomForestSurrogate(xys, zs, lower_bound, upper_bound) +RandomForest = RandomForestSurrogate(xys, zs, lower_bound, upper_bound) ``` ```@example RandomForestSurrogateND p1 = surface(x, y, (x, y) -> RandomForest([x y])) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide p2 = contour(x, y, (x, y) -> RandomForest([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide -plot(p1, p2, title="Surrogate") # hide +scatter!(xs, ys, marker_z = zs) # hide +plot(p1, p2, title = "Surrogate") # hide ``` - ### Optimizing + With our surrogate, we can now search for the minima of the function. Notice how the new sampled points, which were created during the optimization process, are appended to the `xys` array. @@ -125,9 +134,12 @@ This is why its size changes. ```@example RandomForestSurrogateND size(xys) ``` + ```@example RandomForestSurrogateND -surrogate_optimize(bukin6, SRBF(), lower_bound, upper_bound, RandomForest, SobolSample(), maxiters=20) +surrogate_optimize( + bukin6, SRBF(), lower_bound, upper_bound, RandomForest, SobolSample(), maxiters = 20) ``` + ```@example RandomForestSurrogateND size(xys) ``` @@ -137,8 +149,8 @@ p1 = surface(x, y, (x, y) -> RandomForest([x y])) # hide xs = [xy[1] for xy in xys] # hide ys = [xy[2] for xy in xys] # hide zs = bukin6.(xys) # hide -scatter!(xs, ys, zs, marker_z=zs) # hide +scatter!(xs, ys, zs, marker_z = zs) # hide p2 = contour(x, y, (x, y) -> RandomForest([x y])) # hide -scatter!(xs, ys, marker_z=zs) # hide +scatter!(xs, ys, marker_z = zs) # hide plot(p1, p2) # hide ``` diff --git a/docs/src/rosenbrock.md b/docs/src/rosenbrock.md index a882b438b..078612b29 100644 --- a/docs/src/rosenbrock.md +++ b/docs/src/rosenbrock.md @@ -1,10 +1,12 @@ # Rosenbrock function + The Rosenbrock function is defined as: ``f(x) = \sum_{i=1}^{d-1}[ (x_{i+1}-x_i)^2 + (x_i - 1)^2]`` I will treat the 2D version, which is commonly defined as: ``f(x,y) = (1-x)^2 + 100(y-x^2)^2`` Let's import Surrogates and Plots: + ```@example rosen using Surrogates using SurrogatesPolyChaos @@ -13,59 +15,63 @@ default() ``` Define the objective function: + ```@example rosen function f(x) x1 = x[1] x2 = x[2] - return (1-x1)^2 + 100*(x2-x1^2)^2 + return (1 - x1)^2 + 100 * (x2 - x1^2)^2 end ``` Let's plot it: + ```@example rosen n = 100 -lb = [0.0,0.0] -ub = [8.0,8.0] -xys = sample(n,lb,ub,SobolSample()); +lb = [0.0, 0.0] +ub = [8.0, 8.0] +xys = sample(n, lb, ub, SobolSample()); zs = f.(xys); x, y = 0:8, 0:8 -p1 = surface(x, y, (x1,x2) -> f((x1,x2))) +p1 = surface(x, y, (x1, x2) -> f((x1, x2))) xs = [xy[1] for xy in xys] ys = [xy[2] for xy in xys] scatter!(xs, ys, zs) # hide -p2 = contour(x, y, (x1,x2) -> f((x1,x2))) +p2 = contour(x, y, (x1, x2) -> f((x1, x2))) scatter!(xs, ys) -plot(p1, p2, title="True function") +plot(p1, p2, title = "True function") ``` Fitting different surrogates: + ```@example rosen -mypoly = PolynomialChaosSurrogate(xys, zs, lb, ub) +mypoly = PolynomialChaosSurrogate(xys, zs, lb, ub) loba = LobachevskySurrogate(xys, zs, lb, ub) -inver = InverseDistanceSurrogate(xys, zs, lb, ub) +inver = InverseDistanceSurrogate(xys, zs, lb, ub) ``` Plotting: + ```@example rosen p1 = surface(x, y, (x, y) -> mypoly([x y])) -scatter!(xs, ys, zs, marker_z=zs) +scatter!(xs, ys, zs, marker_z = zs) p2 = contour(x, y, (x, y) -> mypoly([x y])) -scatter!(xs, ys, marker_z=zs) -plot(p1, p2, title="Polynomial expansion") +scatter!(xs, ys, marker_z = zs) +plot(p1, p2, title = "Polynomial expansion") ``` ```@example rosen p1 = surface(x, y, (x, y) -> loba([x y])) -scatter!(xs, ys, zs, marker_z=zs) +scatter!(xs, ys, zs, marker_z = zs) p2 = contour(x, y, (x, y) -> loba([x y])) -scatter!(xs, ys, marker_z=zs) -plot(p1, p2, title="Lobachevsky") +scatter!(xs, ys, marker_z = zs) +plot(p1, p2, title = "Lobachevsky") ``` ```@example rosen p1 = surface(x, y, (x, y) -> inver([x y])) -scatter!(xs, ys, zs, marker_z=zs) +scatter!(xs, ys, zs, marker_z = zs) p2 = contour(x, y, (x, y) -> inver([x y])) -scatter!(xs, ys, marker_z=zs) -plot(p1, p2, title="Inverse distance") +scatter!(xs, ys, marker_z = zs) +plot(p1, p2, title = "Inverse distance") ``` diff --git a/docs/src/samples.md b/docs/src/samples.md index 2a92a9d89..17a25b6d7 100644 --- a/docs/src/samples.md +++ b/docs/src/samples.md @@ -3,39 +3,47 @@ Sampling methods are provided by the [QuasiMonteCarlo package](https://docs.sciml.ai/QuasiMonteCarlo/stable/). The syntax for sampling in an interval or region is the following: + ``` sample(n,lb,ub,S::SamplingAlgorithm) ``` + where lb and ub are, respectively, the lower and upper bounds. There are many sampling algorithms to choose from: -* Grid sample + - Grid sample + ``` GridSample{T} sample(n,lb,ub,S::GridSample) ``` -* Uniform sample + - Uniform sample + ``` sample(n,lb,ub,::RandomSample) ``` -* Sobol sample + - Sobol sample + ``` sample(n,lb,ub,::SobolSample) ``` -* Latin Hypercube sample + - Latin Hypercube sample + ``` sample(n,lb,ub,::LatinHypercubeSample) ``` -* Low Discrepancy sample + - Low Discrepancy sample + ``` sample(n,lb,ub,S::HaltonSample) ``` -* Sample on section + - Sample on section + ``` SectionSample sample(n,lb,ub,S::SectionSample) @@ -45,10 +53,11 @@ sample(n,lb,ub,S::SectionSample) Adding a new sampling method is a two- step process: -1. Adding a new SamplingAlgorithm type -2. Overloading the sample function with the new type. + 1. Adding a new SamplingAlgorithm type + 2. Overloading the sample function with the new type. **Example** + ``` struct NewAmazingSamplingAlgorithm{OPTIONAL} <: QuasiMonteCarlo.SamplingAlgorithm end diff --git a/docs/src/secondorderpoly.md b/docs/src/secondorderpoly.md index 790f34d30..eceb84ac5 100644 --- a/docs/src/secondorderpoly.md +++ b/docs/src/secondorderpoly.md @@ -14,30 +14,32 @@ default() ## Sampling ```@example second_order_tut -f = x -> 3*sin(x) + 10/x +f = x -> 3 * sin(x) + 10 / x lb = 3.0 ub = 6.0 n = 10 -x = sample(n,lb,ub,HaltonSample()) +x = sample(n, lb, ub, HaltonSample()) y = f.(x) -scatter(x, y, label="Sampled points", xlims=(lb, ub)) -plot!(f, label="True function", xlims=(lb, ub)) +scatter(x, y, label = "Sampled points", xlims = (lb, ub)) +plot!(f, label = "True function", xlims = (lb, ub)) ``` ## Building the surrogate + ```@example second_order_tut sec = SecondOrderPolynomialSurrogate(x, y, lb, ub) -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lb, ub)) -plot!(f, label="True function", xlims=(lb, ub)) -plot!(sec, label="Surrogate function", xlims=(lb, ub)) +plot(x, y, seriestype = :scatter, label = "Sampled points", xlims = (lb, ub)) +plot!(f, label = "True function", xlims = (lb, ub)) +plot!(sec, label = "Surrogate function", xlims = (lb, ub)) ``` ## Optimizing ```@example second_order_tut @show surrogate_optimize(f, SRBF(), lb, ub, sec, SobolSample()) -scatter(x, y, label="Sampled points") -plot!(f, label="True function", xlims=(lb, ub)) -plot!(sec, label="Surrogate function", xlims=(lb, ub)) +scatter(x, y, label = "Sampled points") +plot!(f, label = "True function", xlims = (lb, ub)) +plot!(sec, label = "Surrogate function", xlims = (lb, ub)) ``` + The optimization method successfully found the minimum. diff --git a/docs/src/sphere_function.md b/docs/src/sphere_function.md index 7c6257d43..dce028f93 100644 --- a/docs/src/sphere_function.md +++ b/docs/src/sphere_function.md @@ -5,6 +5,7 @@ The sphere function of dimension d is defined as: with lower bound -10 and upper bound 10. Let's import Surrogates and Plots: + ```@example sphere_function using Surrogates using Plots @@ -12,44 +13,52 @@ default() ``` Define the objective function: + ```@example sphere_function function sphere_function(x) - return sum(x.^2) + return sum(x .^ 2) end ``` The 1D case is just a simple parabola, let's plot it: + ```@example sphere_function n = 20 lb = -10 ub = 10 -x = sample(n,lb,ub,SobolSample()) +x = sample(n, lb, ub, SobolSample()) y = sphere_function.(x) xs = lb:0.001:ub -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lb, ub), ylims=(-2, 120), legend=:top) -plot!(xs,sphere_function.(xs), label="True function", legend=:top) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lb, ub), ylims = (-2, 120), legend = :top) +plot!(xs, sphere_function.(xs), label = "True function", legend = :top) ``` Fitting RadialSurrogate with different radial basis: + ```@example sphere_function -rad_1d_linear = RadialBasis(x,y,lb,ub) -rad_1d_cubic = RadialBasis(x,y,lb,ub,rad = cubicRadial()) -rad_1d_multiquadric = RadialBasis(x,y,lb,ub, rad = multiquadricRadial()) -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lb, ub), ylims=(-2, 120), legend=:top) -plot!(xs,sphere_function.(xs), label="True function", legend=:top) -plot!(xs, rad_1d_linear.(xs), label="Radial surrogate with linear", legend=:top) -plot!(xs, rad_1d_cubic.(xs), label="Radial surrogate with cubic", legend=:top) -plot!(xs, rad_1d_multiquadric.(xs), label="Radial surrogate with multiquadric", legend=:top) +rad_1d_linear = RadialBasis(x, y, lb, ub) +rad_1d_cubic = RadialBasis(x, y, lb, ub, rad = cubicRadial()) +rad_1d_multiquadric = RadialBasis(x, y, lb, ub, rad = multiquadricRadial()) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lb, ub), ylims = (-2, 120), legend = :top) +plot!(xs, sphere_function.(xs), label = "True function", legend = :top) +plot!(xs, rad_1d_linear.(xs), label = "Radial surrogate with linear", legend = :top) +plot!(xs, rad_1d_cubic.(xs), label = "Radial surrogate with cubic", legend = :top) +plot!(xs, rad_1d_multiquadric.(xs), + label = "Radial surrogate with multiquadric", legend = :top) ``` Fitting Lobachevsky Surrogate with different values of hyperparameter alpha: + ```@example sphere_function -loba_1 = LobachevskySurrogate(x,y,lb,ub) -loba_2 = LobachevskySurrogate(x,y,lb,ub,alpha = 1.5, n = 6) -loba_3 = LobachevskySurrogate(x,y,lb,ub,alpha = 0.3, n = 6) -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lb, ub), ylims=(-2, 120), legend=:top) -plot!(xs,sphere_function.(xs), label="True function", legend=:top) -plot!(xs, loba_1.(xs), label="Lobachevsky surrogate 1", legend=:top) -plot!(xs, loba_2.(xs), label="Lobachevsky surrogate 2", legend=:top) -plot!(xs, loba_3.(xs), label="Lobachevsky surrogate 3", legend=:top) +loba_1 = LobachevskySurrogate(x, y, lb, ub) +loba_2 = LobachevskySurrogate(x, y, lb, ub, alpha = 1.5, n = 6) +loba_3 = LobachevskySurrogate(x, y, lb, ub, alpha = 0.3, n = 6) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lb, ub), ylims = (-2, 120), legend = :top) +plot!(xs, sphere_function.(xs), label = "True function", legend = :top) +plot!(xs, loba_1.(xs), label = "Lobachevsky surrogate 1", legend = :top) +plot!(xs, loba_2.(xs), label = "Lobachevsky surrogate 2", legend = :top) +plot!(xs, loba_3.(xs), label = "Lobachevsky surrogate 3", legend = :top) ``` diff --git a/docs/src/surrogate.md b/docs/src/surrogate.md index 0260fc781..7793f24a6 100644 --- a/docs/src/surrogate.md +++ b/docs/src/surrogate.md @@ -1,57 +1,67 @@ # Surrogate + Every surrogate has a different definition depending on the parameters needed. However, they have in common: -1. ```add_point!(::AbstractSurrogate,x_new,y_new)``` -2. ```AbstractSurrogate(value)``` -The first function adds a sample point to the surrogate, thus changing the internal -coefficients. The second one calculates the approximation at value. + 1. `add_point!(::AbstractSurrogate,x_new,y_new)` + 2. `AbstractSurrogate(value)` + The first function adds a sample point to the surrogate, thus changing the internal + coefficients. The second one calculates the approximation at value. + + - Linear surrogate -* Linear surrogate ```@docs LinearSurrogate(x,y,lb,ub) ``` -* Radial basis function surrogate + - Radial basis function surrogate + ```@docs RadialBasis(x, y, lb, ub; rad::RadialFunction = linearRadial, scale_factor::Real=1.0, sparse = false) ``` -* Kriging surrogate + - Kriging surrogate + ```@docs Kriging(x,y,p,theta) ``` -* Lobachevsky surrogate + - Lobachevsky surrogate + ```@docs LobachevskySurrogate(x,y,lb,ub; alpha = collect(one.(x[1])),n::Int = 4, sparse = false) lobachevsky_integral(loba::LobachevskySurrogate,lb,ub) ``` -* Support vector machine surrogate, requires `using LIBSVM` and `using SurrogatesSVM` + - Support vector machine surrogate, requires `using LIBSVM` and `using SurrogatesSVM` + ``` SVMSurrogate(x,y,lb::Number,ub::Number) ``` -* Random forest surrogate, requires `using XGBoost` and `using SurrogatesRandomForest` + - Random forest surrogate, requires `using XGBoost` and `using SurrogatesRandomForest` + ``` RandomForestSurrogate(x,y,lb,ub;num_round::Int = 1) ``` -* Neural network surrogate, requires `using Flux` and `using SurrogatesFlux` + - Neural network surrogate, requires `using Flux` and `using SurrogatesFlux` + ``` NeuralSurrogate(x,y,lb,ub; model = Chain(Dense(length(x[1]),1), first), loss = (x,y) -> Flux.mse(model(x), y),opt = Descent(0.01),n_echos::Int = 1) ``` # Creating another surrogate + It's great that you want to add another surrogate to the library! You will need to: -1. Define a new mutable struct and a constructor function -2. Define add\_point!(your\_surrogate::AbstractSurrogate,x\_new,y\_new) -3. Define your\_surrogate(value) for the approximation + 1. Define a new mutable struct and a constructor function + 2. Define add\_point!(your\_surrogate::AbstractSurrogate,x\_new,y\_new) + 3. Define your\_surrogate(value) for the approximation **Example** + ``` mutable struct NewSurrogate{X,Y,L,U,C,A,B} <: AbstractSurrogate x::X diff --git a/docs/src/tensor_prod.md b/docs/src/tensor_prod.md index 832603cc4..0288e8214 100644 --- a/docs/src/tensor_prod.md +++ b/docs/src/tensor_prod.md @@ -1,8 +1,10 @@ # Tensor product function + The tensor product function is defined as: ``f(x) = \prod_{i=1}^d \cos(a\pi x_i)`` Let's import Surrogates and Plots: + ```@example tensor using Surrogates using Plots @@ -10,10 +12,11 @@ default() ``` Define the 1D objective function: + ```@example tensor function f(x) - a = 0.5; - return cos(a*pi*x) + a = 0.5 + return cos(a * pi * x) end ``` @@ -25,16 +28,18 @@ a = 0.5 x = sample(n, lb, ub, SobolSample()) y = f.(x) xs = lb:0.001:ub -scatter(x, y, label="Sampled points", xlims=(lb, ub), ylims=(-1, 1), legend=:top) -plot!(xs, f.(xs), label="True function", legend=:top) +scatter(x, y, label = "Sampled points", xlims = (lb, ub), ylims = (-1, 1), legend = :top) +plot!(xs, f.(xs), label = "True function", legend = :top) ``` Fitting and plotting different surrogates: + ```@example tensor loba_1 = LobachevskySurrogate(x, y, lb, ub) krig = Kriging(x, y, lb, ub) -scatter(x, y, label="Sampled points", xlims=(lb, ub), ylims=(-2.5, 2.5), legend=:bottom) -plot!(xs,f.(xs), label="True function", legend=:top) -plot!(xs, loba_1.(xs), label="Lobachevsky", legend=:top) -plot!(xs, krig.(xs), label="Kriging", legend=:top) +scatter( + x, y, label = "Sampled points", xlims = (lb, ub), ylims = (-2.5, 2.5), legend = :bottom) +plot!(xs, f.(xs), label = "True function", legend = :top) +plot!(xs, loba_1.(xs), label = "Lobachevsky", legend = :top) +plot!(xs, krig.(xs), label = "Kriging", legend = :top) ``` diff --git a/docs/src/tutorials.md b/docs/src/tutorials.md index 3dd7d8fac..b08a0b2c2 100644 --- a/docs/src/tutorials.md +++ b/docs/src/tutorials.md @@ -1,16 +1,17 @@ ## Surrogates 101 + Let's start with something easy to get our hands dirty. I want to build a surrogate for ``f(x) = \log(x) \cdot x^2+x^3``. Let's choose the radial basis surrogate. ```@example using Surrogates -f = x -> log(x)*x^2+x^3 +f = x -> log(x) * x^2 + x^3 lb = 1.0 ub = 10.0 -x = sample(50,lb,ub,SobolSample()) +x = sample(50, lb, ub, SobolSample()) y = f.(x) -my_radial_basis = RadialBasis(x,y,lb,ub) +my_radial_basis = RadialBasis(x, y, lb, ub) #I want an approximation at 5.4 approx = my_radial_basis(5.4) @@ -21,18 +22,19 @@ Let's now see an example in 2D. ```@example using Surrogates using LinearAlgebra -f = x -> x[1]*x[2] -lb = [1.0,2.0] -ub = [10.0,8.5] -x = sample(50,lb,ub,SobolSample()) +f = x -> x[1] * x[2] +lb = [1.0, 2.0] +ub = [10.0, 8.5] +x = sample(50, lb, ub, SobolSample()) y = f.(x) -my_radial_basis = RadialBasis(x,y,lb,ub) +my_radial_basis = RadialBasis(x, y, lb, ub) #I want an approximation at (1.0,1.4) -approx = my_radial_basis((1.0,1.4)) +approx = my_radial_basis((1.0, 1.4)) ``` ## Kriging standard error + Let's now use the Kriging surrogate, which is a single-output Gaussian process. This surrogate has a nice feature: not only does it approximate the solution at a point, it also calculates the standard error at such a point. @@ -40,30 +42,32 @@ Let's see an example: ```@example kriging using Surrogates -f = x -> exp(x)*x^2+x^3 +f = x -> exp(x) * x^2 + x^3 lb = 0.0 ub = 10.0 -x = sample(50,lb,ub,RandomSample()) +x = sample(50, lb, ub, RandomSample()) y = f.(x) p = 1.9 -my_krig = Kriging(x,y,lb,ub,p=p) +my_krig = Kriging(x, y, lb, ub, p = p) #I want an approximation at 5.4 approx = my_krig(5.4) #I want to find the standard error at 5.4 -std_err = std_error_at_point(my_krig,5.4) +std_err = std_error_at_point(my_krig, 5.4) ``` Let's now optimize the Kriging surrogate using the lower confidence bound method. This is just a one-liner: ```@example kriging -surrogate_optimize(f,LCBS(),lb,ub,my_krig,RandomSample(); maxiters = 10, num_new_samples = 10) +surrogate_optimize( + f, LCBS(), lb, ub, my_krig, RandomSample(); maxiters = 10, num_new_samples = 10) ``` Surrogate optimization methods have two purposes: they both sample the space in unknown regions and look for the minima at the same time. ## Lobachevsky integral + The Lobachevsky surrogate has the nice feature of having a closed formula for its integral, which is something that other surrogates are missing. Let's compare it with QuadGK: @@ -71,25 +75,25 @@ Let's compare it with QuadGK: ```@example using Surrogates using QuadGK -obj = x -> 3*x + log(x) +obj = x -> 3 * x + log(x) a = 1.0 b = 4.0 -x = sample(2000,a,b,SobolSample()) +x = sample(2000, a, b, SobolSample()) y = obj.(x) alpha = 2.0 n = 6 -my_loba = LobachevskySurrogate(x,y,a,b,alpha=alpha,n=n) +my_loba = LobachevskySurrogate(x, y, a, b, alpha = alpha, n = n) #1D integral -int_1D = lobachevsky_integral(my_loba,a,b) -int = quadgk(obj,a,b) -int_val_true = int[1]-int[2] +int_1D = lobachevsky_integral(my_loba, a, b) +int = quadgk(obj, a, b) +int_val_true = int[1] - int[2] println(int_1D) println(int_val_true) ``` - ## Example of NeuralSurrogate + Basic example of fitting a neural network on a simple function of two variables. ```@example @@ -114,7 +118,8 @@ loss(x, y) = Flux.mse(model(x), y) learning_rate = 0.1 optimizer = Descent(learning_rate) # Simple gradient descent. See Flux documentation for other options. n_epochs = 50 -sgt = NeuralSurrogate(x_train, y_train, bounds..., model=model, loss=loss, opt=optimizer, n_echos=n_epochs) +sgt = NeuralSurrogate(x_train, y_train, bounds..., model = model, + loss = loss, opt = optimizer, n_echos = n_epochs) # Testing the new model x_test = sample(30, bounds..., SobolSample()) diff --git a/docs/src/variablefidelity.md b/docs/src/variablefidelity.md index 180475105..5bd3b1103 100644 --- a/docs/src/variablefidelity.md +++ b/docs/src/variablefidelity.md @@ -13,19 +13,22 @@ default() n = 20 lower_bound = 1.0 upper_bound = 6.0 -x = sample(n,lower_bound,upper_bound,SobolSample()) -f = x -> 1/3*x +x = sample(n, lower_bound, upper_bound, SobolSample()) +f = x -> 1 / 3 * x y = f.(x) -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lower_bound, upper_bound), legend = :top) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) ``` ```@example variablefid -varfid = VariableFidelitySurrogate(x,y,lower_bound,upper_bound) +varfid = VariableFidelitySurrogate(x, y, lower_bound, upper_bound) ``` ```@example variablefid -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) -plot!(varfid, label="Surrogate function", xlims=(lower_bound, upper_bound), legend=:top) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lower_bound, upper_bound), legend = :top) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) +plot!( + varfid, label = "Surrogate function", xlims = (lower_bound, upper_bound), legend = :top) ``` diff --git a/docs/src/water_flow.md b/docs/src/water_flow.md index 885650976..515af9a38 100644 --- a/docs/src/water_flow.md +++ b/docs/src/water_flow.md @@ -14,6 +14,7 @@ default() ``` Define the objective function: + ```@example water function f(x) r_w = x[1] @@ -24,32 +25,31 @@ function f(x) H_l = x[6] L = x[7] K_w = x[8] - log_val = log(r/r_w) - return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l)) + log_val = log(r / r_w) + return (2 * pi * T_u * (H_u - H_l)) / + (log_val * (1 + (2 * L * T_u / (log_val * r_w^2 * K_w)) + T_u / T_l)) end ``` - ```@example water n = 180 d = 8 -lb = [0.05,100,63070,990,63.1,700,1120,9855] -ub = [0.15,50000,115600,1110,116,820,1680,12045] -x = sample(n,lb,ub,SobolSample()) +lb = [0.05, 100, 63070, 990, 63.1, 700, 1120, 9855] +ub = [0.15, 50000, 115600, 1110, 116, 820, 1680, 12045] +x = sample(n, lb, ub, SobolSample()) y = f.(x) n_test = 1000 -x_test = sample(n_test,lb,ub,GoldenSample()); +x_test = sample(n_test, lb, ub, GoldenSample()); y_true = f.(x_test); ``` - ```@example water -my_rad = RadialBasis(x,y,lb,ub) +my_rad = RadialBasis(x, y, lb, ub) y_rad = my_rad.(x_test) -my_poly = PolynomialChaosSurrogate(x,y,lb,ub) +my_poly = PolynomialChaosSurrogate(x, y, lb, ub) y_poly = my_poly.(x_test) -mse_rad = norm(y_true - y_rad,2)/n_test -mse_poly = norm(y_true - y_poly, 2)/n_test +mse_rad = norm(y_true - y_rad, 2) / n_test +mse_poly = norm(y_true - y_poly, 2) / n_test println("MSE Radial: $mse_rad") println("MSE Radial: $mse_poly") ``` diff --git a/docs/src/welded_beam.md b/docs/src/welded_beam.md index c207bce43..c9319f02c 100644 --- a/docs/src/welded_beam.md +++ b/docs/src/welded_beam.md @@ -16,39 +16,39 @@ default() ``` Define the objective function: + ```@example welded function f(x) h = x[1] l = x[2] t = x[3] - a = 6000/(sqrt(2)*h*l) - b = (6000*(14+0.5*l)*sqrt(0.25*(l^2+(h+t)^2)))/(2*(0.707*h*l*(l^2/12 + 0.25*(h+t)^2))) - return (sqrt(a^2+b^2 + l*a*b))/(sqrt(0.25*(l^2+(h+t)^2))) + a = 6000 / (sqrt(2) * h * l) + b = (6000 * (14 + 0.5 * l) * sqrt(0.25 * (l^2 + (h + t)^2))) / + (2 * (0.707 * h * l * (l^2 / 12 + 0.25 * (h + t)^2))) + return (sqrt(a^2 + b^2 + l * a * b)) / (sqrt(0.25 * (l^2 + (h + t)^2))) end ``` - ```@example welded n = 300 d = 3 -lb = [0.125,5.0,5.0] -ub = [1.,10.,10.] -x = sample(n,lb,ub,SobolSample()) +lb = [0.125, 5.0, 5.0] +ub = [1.0, 10.0, 10.0] +x = sample(n, lb, ub, SobolSample()) y = f.(x) n_test = 1000 -x_test = sample(n_test,lb,ub,GoldenSample()); +x_test = sample(n_test, lb, ub, GoldenSample()); y_true = f.(x_test); ``` - ```@example welded -my_rad = RadialBasis(x,y,lb,ub) +my_rad = RadialBasis(x, y, lb, ub) y_rad = my_rad.(x_test) -mse_rad = norm(y_true - y_rad,2)/n_test +mse_rad = norm(y_true - y_rad, 2) / n_test println("MSE Radial: $mse_rad") -my_loba = LobachevskySurrogate(x,y,lb,ub) +my_loba = LobachevskySurrogate(x, y, lb, ub) y_loba = my_loba.(x_test) -mse_rad = norm(y_true - y_loba,2)/n_test +mse_rad = norm(y_true - y_loba, 2) / n_test println("MSE Lobachevsky: $mse_rad") ``` diff --git a/docs/src/wendland.md b/docs/src/wendland.md index 88e975211..b97579cdb 100644 --- a/docs/src/wendland.md +++ b/docs/src/wendland.md @@ -15,7 +15,7 @@ n = 40 lower_bound = 0.0 upper_bound = 1.0 f = x -> exp(-x^2) -x = sample(n,lower_bound,upper_bound,SobolSample()) +x = sample(n, lower_bound, upper_bound, SobolSample()) y = f.(x) ``` @@ -29,11 +29,12 @@ Try it yourself with this function! ```@example wendland my_eps = 0.5 -wend = Wendland(x,y,lower_bound,upper_bound,eps=my_eps) +wend = Wendland(x, y, lower_bound, upper_bound, eps = my_eps) ``` ```@example wendland -plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound), legend=:top) -plot!(f, label="True function", xlims=(lower_bound, upper_bound), legend=:top) -plot!(wend, label="Surrogate function", xlims=(lower_bound, upper_bound), legend=:top) +plot(x, y, seriestype = :scatter, label = "Sampled points", + xlims = (lower_bound, upper_bound), legend = :top) +plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top) +plot!(wend, label = "Surrogate function", xlims = (lower_bound, upper_bound), legend = :top) ``` diff --git a/lib/SurrogatesAbstractGPs/LICENSE.md b/lib/SurrogatesAbstractGPs/LICENSE.md index 62199eff7..27dd257c0 100644 --- a/lib/SurrogatesAbstractGPs/LICENSE.md +++ b/lib/SurrogatesAbstractGPs/LICENSE.md @@ -1,7 +1,7 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > Copyright (c) 2019-20: Ludovico Bessi, Julia Computing. -> +> > Permission is hereby granted, free of charge, to any person obtaining > a copy of this software and associated documentation files (the > "Software"), to deal in the Software without restriction, including @@ -9,10 +9,10 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > distribute, sublicense, and/or sell copies of the Software, and to > permit persons to whom the Software is furnished to do so, subject to > the following conditions: -> +> > The above copyright notice and this permission notice shall be > included in all copies or substantial portions of the Software. -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. diff --git a/lib/SurrogatesFlux/LICENSE.md b/lib/SurrogatesFlux/LICENSE.md index 62199eff7..27dd257c0 100644 --- a/lib/SurrogatesFlux/LICENSE.md +++ b/lib/SurrogatesFlux/LICENSE.md @@ -1,7 +1,7 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > Copyright (c) 2019-20: Ludovico Bessi, Julia Computing. -> +> > Permission is hereby granted, free of charge, to any person obtaining > a copy of this software and associated documentation files (the > "Software"), to deal in the Software without restriction, including @@ -9,10 +9,10 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > distribute, sublicense, and/or sell copies of the Software, and to > permit persons to whom the Software is furnished to do so, subject to > the following conditions: -> +> > The above copyright notice and this permission notice shall be > included in all copies or substantial portions of the Software. -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. diff --git a/lib/SurrogatesFlux/src/SurrogatesFlux.jl b/lib/SurrogatesFlux/src/SurrogatesFlux.jl index f97085fb3..fa58a3f0b 100644 --- a/lib/SurrogatesFlux/src/SurrogatesFlux.jl +++ b/lib/SurrogatesFlux/src/SurrogatesFlux.jl @@ -20,10 +20,9 @@ end """ NeuralSurrogate(x,y,lb,ub,model,loss,opt,n_echos) -- model: Flux layers -- loss: loss function -- opt: optimization function - + - model: Flux layers + - loss: loss function + - opt: optimization function """ function NeuralSurrogate(x, y, lb, ub; model = Chain(Dense(length(x[1]), 1), first), loss = (x, y) -> Flux.mse(model(x), y), opt = Descent(0.01), diff --git a/lib/SurrogatesMOE/LICENSE.md b/lib/SurrogatesMOE/LICENSE.md index 62199eff7..27dd257c0 100644 --- a/lib/SurrogatesMOE/LICENSE.md +++ b/lib/SurrogatesMOE/LICENSE.md @@ -1,7 +1,7 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > Copyright (c) 2019-20: Ludovico Bessi, Julia Computing. -> +> > Permission is hereby granted, free of charge, to any person obtaining > a copy of this software and associated documentation files (the > "Software"), to deal in the Software without restriction, including @@ -9,10 +9,10 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > distribute, sublicense, and/or sell copies of the Software, and to > permit persons to whom the Software is furnished to do so, subject to > the following conditions: -> +> > The above copyright notice and this permission notice shall be > included in all copies or substantial portions of the Software. -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. diff --git a/lib/SurrogatesMOE/src/SurrogatesMOE.jl b/lib/SurrogatesMOE/src/SurrogatesMOE.jl index 4ac49ea9a..4b3369528 100644 --- a/lib/SurrogatesMOE/src/SurrogatesMOE.jl +++ b/lib/SurrogatesMOE/src/SurrogatesMOE.jl @@ -1,10 +1,10 @@ module SurrogatesMOE import Surrogates: AbstractSurrogate, linearRadial, cubicRadial, multiquadricRadial, - thinplateRadial, RadialBasisStructure, RadialBasis, - InverseDistanceSurrogate, Kriging, LobachevskyStructure, - LobachevskySurrogate, NeuralStructure, PolyChaosStructure, - LinearSurrogate, add_point! + thinplateRadial, RadialBasisStructure, RadialBasis, + InverseDistanceSurrogate, Kriging, LobachevskyStructure, + LobachevskySurrogate, NeuralStructure, PolyChaosStructure, + LinearSurrogate, add_point! export MOE @@ -30,6 +30,7 @@ end """ MOE(x, y, expert_types; ndim=1, n_clusters=2) + constructor for MOE; takes in x, y and expert types and returns an MOE struct """ function MOE(x, y, expert_types; ndim = 1, n_clusters = 2, quantile = 10) @@ -68,6 +69,7 @@ end """ (moe::MOE)(val::Number) + predictor for 1D inputs """ function (moe::MOE)(val::Number) @@ -88,7 +90,6 @@ end (moe::MOE)(val) predictor for ndimensional inputs - """ function (moe::MOE)(val) val = collect(val) #to handle inputs that may sometimes be tuples @@ -108,9 +109,10 @@ end """ _cluster_predict(gmm:GMM, X::Matrix) + gmm - a trained Gaussian Mixture Model -X - a matrix of points with dimensions equal to the inputs used for the - training of the model +X - a matrix of points with dimensions equal to the inputs used for the +training of the model Return - Clusters to which each of the points belong to (starts at int 1) @@ -147,6 +149,7 @@ end """ _cluster_values(values, cluster_classifier, num_clusters) + values - a concatenation of input and output values cluster_classifier - a vector of integers representing which cluster each data point belongs to num_clusters - number of clusters @@ -154,16 +157,16 @@ num_clusters - number of clusters output clusters - values grouped by clusters -Ex: +## Ex: + vals = [1.0 2.0; 3.0 4.0; 5.0 6.0; 7.0 8.0; 9.0 10.0] cluster_classifier = [1, 2, 2, 2, 1] num_clusters = 2 clusters = _cluster_values(vals, cluster_classifier, num_clusters) @show clusters #prints values below ---- + [[1.0, 2.0], [9.0, 10.0]] [[3.0, 4.0], [5.0, 6.0], [7.0, 8.0]] - """ function _cluster_values(values, cluster_classifier, num_clusters) num = length(cluster_classifier) @@ -202,7 +205,7 @@ end """ _find_upper_lower_bounds(m::Matrix) - returns upper and lower bounds in vector form +returns upper and lower bounds in vector form """ function _find_upper_lower_bounds(X::Matrix) ub = [] @@ -221,7 +224,6 @@ end """ _find_best_model(clustered_values, clustered_test_values) finds best model for each set of clustered values by validating against the clustered_test_values - """ function _find_best_model(clustered_train_values, clustered_test_values, dim, enabled_expert_types) @@ -248,7 +250,8 @@ function _find_best_model(clustered_train_values, clustered_test_values, dim, # call on _surrogate_builder with clustered_train_vals, enabled expert types, lb, ub - surr_vec = _surrogate_builder(enabled_expert_types, length(enabled_expert_types), x_vec, + surr_vec = _surrogate_builder( + enabled_expert_types, length(enabled_expert_types), x_vec, y_vec, lb, ub) # use the models to find best model after validating against test data and return best model @@ -267,6 +270,7 @@ end """ _surrogate_builder(local_kind, k, x, y, lb, ub) + takes in an array of surrogate types, and number of cluster, builds the surrogates and returns an array of surrogate objects """ @@ -384,6 +388,7 @@ end """ _vector_of_tuples_to_matrix(v) + takes in a vector of tuples or vector of vectors and converts it into a matrix """ function _vector_of_tuples_to_matrix(v) diff --git a/lib/SurrogatesMOE/test/runtests.jl b/lib/SurrogatesMOE/test/runtests.jl index 43f9e8d03..6e094cda2 100644 --- a/lib/SurrogatesMOE/test/runtests.jl +++ b/lib/SurrogatesMOE/test/runtests.jl @@ -28,7 +28,7 @@ Random.seed!(StableRNG(SEED), SEED) RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, sparse = false), RadialBasisStructure(radial_function = cubicRadial(), scale_factor = 1.0, - sparse = false), + sparse = false) ] MOE_1D_RAD_RAD = MOE(x, y, expert_types) @@ -40,7 +40,7 @@ Random.seed!(StableRNG(SEED), SEED) # Krig vs MOE KRIG_1D = Kriging(x, y, lb, ub, p = 1.0, theta = 1.0) expert_types = [InverseDistanceStructure(p = 1.0), - KrigingStructure(p = 1.0, theta = 1.0), + KrigingStructure(p = 1.0, theta = 1.0) ] MOE_1D_INV_KRIG = MOE(x, y, expert_types) MOE_at0 = MOE_1D_INV_KRIG(0.0) @@ -83,7 +83,7 @@ end expert_types = [ KrigingStructure(p = [1.0, 1.0], theta = [1.0, 1.0]), RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, - sparse = false), + sparse = false) ] moe_nd_krig_rad = MOE(x, y, expert_types, ndim = 2, quantile = 5) moe_pred_vals = moe_nd_krig_rad.(x_test) @@ -125,7 +125,7 @@ end RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, sparse = false), LinearStructure(), - InverseDistanceStructure(p = 1.0), + InverseDistanceStructure(p = 1.0) ] moe_nd_3_experts = MOE(x, y, expert_types, ndim = 2, n_clusters = 3) moe_pred_vals = moe_nd_3_experts.(x_test) @@ -137,7 +137,7 @@ end n_echos = 1 expert_types = [ NeuralStructure(model = model, loss = loss, opt = opt, n_echos = n_echos), - LinearStructure(), + LinearStructure() ] moe_nn_ln = MOE(x, y, expert_types, ndim = 2) moe_pred_vals = moe_nn_ln.(x_test) @@ -163,7 +163,7 @@ end RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, sparse = false), RadialBasisStructure(radial_function = cubicRadial(), scale_factor = 1.0, - sparse = false), + sparse = false) ] moe = MOE(x, y, expert_types) add_point!(moe, 0.5, 5.0) @@ -188,7 +188,7 @@ end y = discont_NDIM.(x) expert_types = [InverseDistanceStructure(p = 1.0), RadialBasisStructure(radial_function = linearRadial(), scale_factor = 1.0, - sparse = false), + sparse = false) ] moe_nd_inv_rad = MOE(x, y, expert_types, ndim = 2) add_point!(moe_nd_inv_rad, (0.5, 0.5), sum((0.5, 0.5) .^ 2) + 5) diff --git a/lib/SurrogatesPolyChaos/LICENSE.md b/lib/SurrogatesPolyChaos/LICENSE.md index 62199eff7..27dd257c0 100644 --- a/lib/SurrogatesPolyChaos/LICENSE.md +++ b/lib/SurrogatesPolyChaos/LICENSE.md @@ -1,7 +1,7 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > Copyright (c) 2019-20: Ludovico Bessi, Julia Computing. -> +> > Permission is hereby granted, free of charge, to any person obtaining > a copy of this software and associated documentation files (the > "Software"), to deal in the Software without restriction, including @@ -9,10 +9,10 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > distribute, sublicense, and/or sell copies of the Software, and to > permit persons to whom the Software is furnished to do so, subject to > the following conditions: -> +> > The above copyright notice and this permission notice shall be > included in all copies or substantial portions of the Software. -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. diff --git a/lib/SurrogatesRandomForest/LICENSE.md b/lib/SurrogatesRandomForest/LICENSE.md index 62199eff7..27dd257c0 100644 --- a/lib/SurrogatesRandomForest/LICENSE.md +++ b/lib/SurrogatesRandomForest/LICENSE.md @@ -1,7 +1,7 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > Copyright (c) 2019-20: Ludovico Bessi, Julia Computing. -> +> > Permission is hereby granted, free of charge, to any person obtaining > a copy of this software and associated documentation files (the > "Software"), to deal in the Software without restriction, including @@ -9,10 +9,10 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > distribute, sublicense, and/or sell copies of the Software, and to > permit persons to whom the Software is furnished to do so, subject to > the following conditions: -> +> > The above copyright notice and this permission notice shall be > included in all copies or substantial portions of the Software. -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. diff --git a/lib/SurrogatesRandomForest/src/SurrogatesRandomForest.jl b/lib/SurrogatesRandomForest/src/SurrogatesRandomForest.jl index b3db175d4..498688748 100644 --- a/lib/SurrogatesRandomForest/src/SurrogatesRandomForest.jl +++ b/lib/SurrogatesRandomForest/src/SurrogatesRandomForest.jl @@ -28,7 +28,6 @@ end RandomForestSurrogate(x,y,lb,ub,num_round) Build Random forest surrogate. num_round is the number of trees. - """ function RandomForestSurrogate(x, y, lb, ub; num_round::Int = 1) X = Array{Float64, 2}(undef, length(x), length(x[1])) diff --git a/lib/SurrogatesSVM/LICENSE.md b/lib/SurrogatesSVM/LICENSE.md index 62199eff7..27dd257c0 100644 --- a/lib/SurrogatesSVM/LICENSE.md +++ b/lib/SurrogatesSVM/LICENSE.md @@ -1,7 +1,7 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > Copyright (c) 2019-20: Ludovico Bessi, Julia Computing. -> +> > Permission is hereby granted, free of charge, to any person obtaining > a copy of this software and associated documentation files (the > "Software"), to deal in the Software without restriction, including @@ -9,10 +9,10 @@ The Surrogates.jl package is licensed under the MIT "Expat" License: > distribute, sublicense, and/or sell copies of the Software, and to > permit persons to whom the Software is furnished to do so, subject to > the following conditions: -> +> > The above copyright notice and this permission notice shall be > included in all copies or substantial portions of the Software. -> +> > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. diff --git a/src/Earth.jl b/src/Earth.jl index fb620b745..d3fab9769 100644 --- a/src/Earth.jl +++ b/src/Earth.jl @@ -229,7 +229,7 @@ function _forward_pass_nd(x, y, n_max_terms, rel_res_error, maxiters) for a in 1:n val_a = sum(coeff[b] * prod([new_basis[b][c](x[a][c]) for c in 1:d]) - for b in 1:bas_len) + intercept + for b in 1:bas_len) + intercept new_sse = new_sse + (y[a] - val_a)^2 end end @@ -271,7 +271,7 @@ function _backward_pass_nd(x, y, n_min_terms, basis, penalty, rel_GCV) sse = zero(y[1]) for a in 1:n val_a = sum(coeff[b] * prod([basis[b][c](x[a][c]) for c in 1:d]) - for b in 1:base_len) + + for b in 1:base_len) + intercept sse = sse + (y[a] - val_a)^2 end @@ -295,7 +295,7 @@ function _backward_pass_nd(x, y, n_min_terms, basis, penalty, rel_GCV) current_base_len = num_terms - 1 for a in 1:n val_a = sum(coeff[b] * prod([current_basis[b][c](x[a][c]) for c in 1:d]) - for b in 1:current_base_len) + intercept + for b in 1:current_base_len) + intercept new_sse = new_sse + (y[a] - val_a)^2 end curr_effect_num_params = current_base_len + penalty * (current_base_len - 1) / 2 diff --git a/src/GEK.jl b/src/GEK.jl index 75da4ca3b..790b2ba0a 100644 --- a/src/GEK.jl +++ b/src/GEK.jl @@ -202,7 +202,7 @@ function (k::GEK)(val) return k.mu + sum(k.b[i] * exp(-sum(k.theta[j] * norm(val[j] - k.x[i][j])^k.p[j] for j in 1:d)) - for i in 1:n) + for i in 1:n) end function GEK(x, y, lb, ub; p = collect(one.(x[1])), theta = collect(one.(x[1]))) diff --git a/src/GEKPLS.jl b/src/GEKPLS.jl index 1e26bf26f..d8716d52d 100644 --- a/src/GEKPLS.jl +++ b/src/GEKPLS.jl @@ -3,8 +3,8 @@ using Statistics """ GEKPLS(x, y, x_matrix, y_matrix, grads, xlimits, delta_x, extra_points, n_comp, beta, gamma, theta, - reduced_likelihood_function_value, - X_offset, X_scale, X_after_std, pls_mean_reshaped, y_mean, y_std) +reduced_likelihood_function_value, +X_offset, X_scale, X_after_std, pls_mean_reshaped, y_mean, y_std) """ mutable struct GEKPLS{T, X, Y} <: AbstractSurrogate x::X @@ -30,6 +30,7 @@ end """ bounds_error(x, xl) + Checks if input x values are within range of upper and lower bounds """ function bounds_error(x, xl) @@ -49,15 +50,16 @@ end GEKPLS(X, y, grads, n_comp, delta_x, lb, ub, extra_points, theta) Constructor for GEKPLS Struct -- x_vec: vector of tuples with x values -- y_vec: vector of floats with outputs -- grads_vec: gradients associated with each of the X points -- n_comp: number of components -- lb: lower bounds -- ub: upper bounds -- delta_x: step size while doing Taylor approximation -- extra_points: number of points to consider -- theta: initial expected variance of PLS regression components + + - x_vec: vector of tuples with x values + - y_vec: vector of floats with outputs + - grads_vec: gradients associated with each of the X points + - n_comp: number of components + - lb: lower bounds + - ub: upper bounds + - delta_x: step size while doing Taylor approximation + - extra_points: number of points to consider + - theta: initial expected variance of PLS regression components """ function GEKPLS(x_vec, y_vec, grads_vec, n_comp, delta_x, lb, ub, extra_points, theta) xlimits = hcat(lb, ub) @@ -73,7 +75,8 @@ function GEKPLS(x_vec, y_vec, grads_vec, n_comp, delta_x, lb, ub, extra_points, pls_mean, X_after_PLS, y_after_PLS = _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) - X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = standardization(X_after_PLS, + X_after_std, y_after_std, X_offset, y_mean, X_scale, y_std = standardization( + X_after_PLS, y_after_PLS) D, ij = cross_distances(X_after_std) pls_mean_reshaped = reshape(pls_mean, (size(X, 2), n_comp)) @@ -91,7 +94,7 @@ end """ (g::GEKPLS)(X_test) - + Take in a set of one or more points and provide predicted approximate outputs (predictor). """ function (g::GEKPLS)(x_vec) @@ -137,13 +140,15 @@ function add_point!(g::GEKPLS, x_tup, y_val, grad_tup) g.num_components, g.grads, g.delta, g.xl, g.extra_points) - g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = standardization(X_after_PLS, + g.X_after_std, y_after_std, g.X_offset, g.y_mean, g.X_scale, g.y_std = standardization( + X_after_PLS, y_after_PLS) D, ij = cross_distances(g.X_after_std) g.pls_mean = reshape(pls_mean, (size(g.x_matrix, 2), g.num_components)) d = componentwise_distance_PLS(D, "squar_exp", g.num_components, g.pls_mean) nt, nd = size(X_after_PLS) - g.beta, g.gamma, g.reduced_likelihood_function_value = _reduced_likelihood_function(g.theta, + g.beta, g.gamma, g.reduced_likelihood_function_value = _reduced_likelihood_function( + g.theta, "squar_exp", d, nt, @@ -153,21 +158,25 @@ end """ _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) -Gradient-enhanced PLS-coefficients. + +## Gradient-enhanced PLS-coefficients. + Parameters ----------- -- X: [n_obs,dim] - The input variables. -- y: [n_obs,ny] - The output variable -- n_comp: int - Number of principal components used. -- gradients: - The gradient values. Matrix size (n_obs,dim) -- delta_x: real - The step used in the First Order Taylor Approximation -- xlimits: [dim, 2]- The upper and lower var bounds. -- extra_points: int - The number of extra points per each training point. -Returns -------- -- Coeff_pls: [dim, n_comp] - The PLS-coefficients. -- X: Concatenation of XX: [extra_points*nt, dim] - Extra points added (when extra_points > 0) and X -- y: Concatenation of yy[extra_points*nt, 1]- Extra points added (when extra_points > 0) and y + + - X: [n_obs,dim] - The input variables. + - y: [n_obs,ny] - The output variable + - n_comp: int - Number of principal components used. + - gradients: - The gradient values. Matrix size (n_obs,dim) + - delta_x: real - The step used in the First Order Taylor Approximation + - xlimits: [dim, 2]- The upper and lower var bounds. + - extra_points: int - The number of extra points per each training point. + Returns + +* * * + + - Coeff_pls: [dim, n_comp] - The PLS-coefficients. + - X: Concatenation of XX: [extra_points*nt, dim] - Extra points added (when extra_points > 0) and X + - y: Concatenation of yy[extra_points*nt, 1]- Extra points added (when extra_points > 0) and y """ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) @@ -185,14 +194,14 @@ function _ge_compute_pls(X, y, n_comp, grads, delta_x, xlimits, extra_points) bb_vals = circshift(boxbehnken(dim, 1), 1) else bb_vals = [0.0 0.0; #center - 1.0 0.0; #right - 0.0 1.0; #up - -1.0 0.0; #left - 0.0 -1.0; #down - 1.0 1.0; #right up - -1.0 1.0; #left up - -1.0 -1.0; #left down - 1.0 -1.0] + 1.0 0.0; #right + 0.0 1.0; #up + -1.0 0.0; #left + 0.0 -1.0; #down + 1.0 1.0; #right up + -1.0 1.0; #left up + -1.0 -1.0; #left down + 1.0 -1.0] end _X = zeros((size(bb_vals)[1], dim)) _y = zeros((size(bb_vals)[1], 1)) @@ -307,20 +316,18 @@ end We subtract the mean from each variable. Then, we divide the values of each variable by its standard deviation. -Parameters ----------- +## Parameters X - The input variables. y - The output variable. -Returns -------- +## Returns X: [n_obs, dim] - The standardized input matrix. +The standardized input matrix. y: [n_obs, 1] - The standardized output vector. +The standardized output vector. X_offset: The mean (or the min if scale_X_to_unit=True) of each input variable. @@ -329,7 +336,6 @@ y_mean: The mean of the output variable. X_scale: The standard deviation of each input variable. y_std: The standard deviation of the output variable. - """ function standardization(X, y) #Equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L21 @@ -348,19 +354,20 @@ end Computes the nonzero componentwise cross-distances between the vectors in X -Parameters ----------- +## Parameters X: [n_obs, dim] -Returns -------- +## Returns + D: [n_obs * (n_obs - 1) / 2, dim] - - The cross-distances between the vectors in X. + + - The cross-distances between the vectors in X. ij: [n_obs * (n_obs - 1) / 2, 2] - - The indices i and j of the vectors in X associated to the cross- - distances in D. + + - The indices i and j of the vectors in X associated to the cross- + distances in D. """ function cross_distances(X) # equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L86 @@ -409,7 +416,6 @@ end D_corr: [n_obs * (n_obs - 1) / 2, n_comp] - The componentwise cross-spatial-correlation-distance between the vectors in X. - """ function componentwise_distance_PLS(D, corr, n_comp, coeff_pls) @@ -430,17 +436,17 @@ function componentwise_distance_PLS(D, corr, n_comp, coeff_pls) end """ -Squared exponential correlation model. +## Squared exponential correlation model. + Equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L604 Parameters: ------------ + theta : Hyperparameters of the correlation model d: componentwise distances from componentwise_distance_PLS -Returns: --------- -r: array containing the values of the autocorrelation model +## Returns: +r: array containing the values of the autocorrelation model """ function squar_exp(theta, d) n_components = size(d)[2] @@ -450,6 +456,7 @@ end """ differences(X, Y) + return differences between two arrays given an input like this: @@ -461,12 +468,12 @@ diff = differences(X,Y) We get an output (diff) that looks like this: [ 0. -1. -2. - -3. -4. -5. - -6. -7. -8. - 1. 0. -1. - -2. -3. -4. - -5. -6. -7.] +-3. -4. -5. +-6. -7. -8. + 1. 0. -1. + -2. -3. -4. + -5. -6. -7.] """ function differences(X, Y) #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/utils/kriging_utils.py#L392 @@ -478,14 +485,14 @@ end """ _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, noise = 0.0) - + Compute the reduced likelihood function value and other coefficients necessary for prediction This function is a loose translation of SMT code from https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/surrogate_models/krg_based.py#L247 It determines the BLUP parameters and evaluates the reduced likelihood function for the given theta. -Parameters ----------- +## Parameters + theta: array containing the parameters at which the Gaussian Process model parameters should be determined. kernel_type: name of the correlation function. d: The componentwise cross-spatial-correlation-distance between the vectors in X. @@ -494,14 +501,13 @@ ij: The indices i and j of the vectors in X associated to the cross-distances in y_norma: Standardized y values noise: noise hyperparameter - increasing noise reduces reduced_likelihood_function_value +## Returns -Returns -------- reduced_likelihood_function_value: real - - The value of the reduced likelihood function associated with the given autocorrelation parameters theta. -beta: Generalized least-squares regression weights -gamma: Gaussian Process weights. + - The value of the reduced likelihood function associated with the given autocorrelation parameters theta. + beta: Generalized least-squares regression weights + gamma: Gaussian Process weights. """ function _reduced_likelihood_function(theta, kernel_type, d, nt, ij, y_norma, noise = 0.0) #equivalent of https://github.com/SMTorg/smt/blob/4a4df255b9259965439120091007f9852f41523e/smt/surrogate_models/krg_based.py#L247 diff --git a/src/InverseDistanceSurrogate.jl b/src/InverseDistanceSurrogate.jl index 856926bd0..d241e0565 100644 --- a/src/InverseDistanceSurrogate.jl +++ b/src/InverseDistanceSurrogate.jl @@ -5,7 +5,6 @@ The inverse distance weighting model is an interpolating method and the unknown points are calculated with a weighted average of the sampling points. p is a positive real number called the power parameter. p > 1 is needed for the derivative to be continuous. - """ mutable struct InverseDistanceSurrogate{X, Y, L, U, P} <: AbstractSurrogate x::X @@ -28,9 +27,9 @@ function (inverSurr::InverseDistanceSurrogate)(val) else if length(inverSurr.lb) == 1 num = sum(inverSurr.y[i] * (norm(val .- inverSurr.x[i]))^(-inverSurr.p) - for i in 1:length(inverSurr.x)) + for i in 1:length(inverSurr.x)) den = sum(norm(val .- inverSurr.x[i])^(-inverSurr.p) - for i in 1:length(inverSurr.x)) + for i in 1:length(inverSurr.x)) return num / den else βᵢ = [norm(val .- inverSurr.x[i])^(-inverSurr.p) for i in 1:length(inverSurr.x)] diff --git a/src/Kriging.jl b/src/Kriging.jl index f48306069..4fd0be3e2 100644 --- a/src/Kriging.jl +++ b/src/Kriging.jl @@ -31,7 +31,7 @@ function (k::Kriging)(val) return k.mu + sum(k.b[i] * exp(-sum(k.theta[j] * norm(val[j] - k.x[i][j])^k.p[j] for j in 1:d)) - for i in 1:n) + for i in 1:n) end """ @@ -46,12 +46,12 @@ function std_error_at_point(k::Kriging, val) d = length(k.x[1]) r = zeros(eltype(k.x[1]), n, 1) r = [let - sum = zero(eltype(k.x[1])) - for l in 1:d - sum = sum + k.theta[l] * norm(val[l] - k.x[i][l])^(k.p[l]) - end - exp(-sum) - end + sum = zero(eltype(k.x[1])) + for l in 1:d + sum = sum + k.theta[l] * norm(val[l] - k.x[i][l])^(k.p[l]) + end + exp(-sum) + end for i in 1:n] one = ones(eltype(k.x[1]), n, 1) @@ -98,8 +98,9 @@ Constructor for type Kriging. #Arguments: -(x,y): sampled points -p: value between 0 and 2 modelling the - smoothness of the function being approximated, 0-> rough 2-> C^infinity -- theta: value > 0 modeling how much the function is changing in the i-th variable. +smoothness of the function being approximated, 0-> rough 2-> C^infinity + + - theta: value > 0 modeling how much the function is changing in the i-th variable. """ function Kriging(x, y, lb::Number, ub::Number; p = 2.0, theta = 0.5 / max(1e-6 * abs(ub - lb), std(x))^p) @@ -160,12 +161,12 @@ end Constructor for Kriging surrogate. -- (x,y): sampled points -- p: array of values 0<=p<2 modeling the - smoothness of the function being approximated in the i-th variable. - low p -> rough, high p -> smooth -- theta: array of values > 0 modeling how much the function is - changing in the i-th variable. + - (x,y): sampled points + - p: array of values 0<=p<2 modeling the + smoothness of the function being approximated in the i-th variable. + low p -> rough, high p -> smooth + - theta: array of values > 0 modeling how much the function is + changing in the i-th variable. """ function Kriging(x, y, lb, ub; p = 2.0 .* collect(one.(x[1])), theta = [0.5 / max(1e-6 * norm(ub .- lb), std(x_i[i] for x_i in x))^p[i] @@ -194,12 +195,12 @@ function _calc_kriging_coeffs(x, y, p, theta) d = length(x[1]) R = [let - sum = zero(eltype(x[1])) - for l in 1:d - sum = sum + theta[l] * norm(x[i][l] - x[j][l])^p[l] - end - exp(-sum) - end + sum = zero(eltype(x[1])) + for l in 1:d + sum = sum + theta[l] * norm(x[i][l] - x[j][l])^p[l] + end + exp(-sum) + end for j in 1:n, i in 1:n] # Estimate nugget based on maximum allowed condition number @@ -243,7 +244,6 @@ end Adds the new point and its respective value to the sample points. Warning: If you are just adding a single point, you have to wrap it with []. Returns the updated Kriging model. - """ function add_point!(k::Kriging, new_x, new_y) if new_x in k.x diff --git a/src/LinearSurrogate.jl b/src/LinearSurrogate.jl index 53be7ba88..496bc319c 100644 --- a/src/LinearSurrogate.jl +++ b/src/LinearSurrogate.jl @@ -60,7 +60,6 @@ end LinearSurrogate(x,y,lb,ub) Builds a linear surrogate using GLM.jl - """ function LinearSurrogate(x, y, lb, ub) #X = Array{eltype(x[1]),2}(undef,length(x),length(x[1])) diff --git a/src/Lobachevsky.jl b/src/Lobachevsky.jl index 30c890afb..34b532253 100644 --- a/src/Lobachevsky.jl +++ b/src/Lobachevsky.jl @@ -44,7 +44,8 @@ end """ Lobachevsky interpolation, suggested parameters: 0 <= alpha <= 4, n must be even. """ -function LobachevskySurrogate(x, y, lb::Number, ub::Number; alpha::Number = 1.0, n::Int = 4, +function LobachevskySurrogate( + x, y, lb::Number, ub::Number; alpha::Number = 1.0, n::Int = 4, sparse = false) if alpha > 4 || alpha < 0 error("Alpha must be between 0 and 4") @@ -61,7 +62,7 @@ function (loba::LobachevskySurrogate)(val::Number) _check_dimension(loba, val) return sum(loba.coeff[j] * phi_nj1D(val, loba.x[j], loba.alpha, loba.n) - for j in 1:length(loba.x)) + for j in 1:length(loba.x)) end function phi_njND(point, x, alpha, n) @@ -101,7 +102,7 @@ function (loba::LobachevskySurrogate)(val) # Check to make sure dimensions of input matches expected dimension of surrogate _check_dimension(loba, val) return sum(loba.coeff[j] * phi_njND(val, loba.x[j], loba.alpha, loba.n) - for j in 1:length(loba.x)) + for j in 1:length(loba.x)) end function add_point!(loba::LobachevskySurrogate, x_new, y_new) @@ -147,7 +148,6 @@ end lobachevsky_integral(loba::LobachevskySurrogate,lb,ub) Calculates the integral of the Lobachevsky surrogate, which has a closed form. - """ function lobachevsky_integral(loba::LobachevskySurrogate, lb, ub) d = length(lb) @@ -168,7 +168,6 @@ end lobachevsky_integrate_dimension(loba::LobachevskySurrogate,lb,ub,dimension) Integrating the surrogate on selected dimension dim - """ function lobachevsky_integrate_dimension(loba::LobachevskySurrogate, lb, ub, dim::Int) gamma_d = zero(loba.coeff[1]) diff --git a/src/Optimization.jl b/src/Optimization.jl index f42ba2d73..7f5fa0322 100755 --- a/src/Optimization.jl +++ b/src/Optimization.jl @@ -76,7 +76,7 @@ the weighted-distance merit function: ``merit(x) = ws(x) + (1-w)(1-d(x))`` -where `` 0 \\leq w \\leq 1 ``. +where ``0 \\leq w \\leq 1``. That is, we want a small function value prediction and a large minimum distance from the previously evaluated points. The weight w is commonly cycled between @@ -848,8 +848,6 @@ Under a Gaussian process (GP) prior, the goal is to maximize expected improvement: ``EI(x) := E[max(f_{best}-f(x),0)`` - - """ function surrogate_optimize(obj::Function, ::EI, lb, ub, krig, sample_type::SamplingAlgorithm; maxiters = 100, @@ -992,7 +990,6 @@ surrogate_optimize(obj::Function,::DYCORS,lb::Number,ub::Number,surr1::AbstractS DYCORS optimization method in 1D, following closely: Combining radial basis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization". - """ function surrogate_optimize(obj::Function, ::DYCORS, lb::Number, ub::Number, surr1::AbstractSurrogate, sample_type::SamplingAlgorithm; @@ -1110,7 +1107,6 @@ end """ surrogate_optimize(obj::Function,::DYCORS,lb::Number,ub::Number,surr1::AbstractSurrogate,sample_type::SamplingAlgorithm;maxiters=100,num_new_samples=100) - This is an implementation of the DYCORS strategy by Regis and Shoemaker: Rommel G Regis and Christine A Shoemaker. Combining radial basis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization. @@ -1307,6 +1303,7 @@ SOP Surrogate optimization method, following closely the following papers: - SOP: parallel surrogate global optimization with Pareto center selection for computationally expensive single objective problems by Tipaluck Krityakierne - Multiobjective Optimization Using Evolutionary Algorithms by Kalyan Deb + #Suggested number of new_samples = min(500*d,5000) """ function surrogate_optimize(obj::Function, sop1::SOP, lb::Number, ub::Number, @@ -1739,8 +1736,10 @@ function _nonDominatedSorting(arr::Array{Float64, 2}) ind::Array{Int64, 1} = collect(1:size(arr, 1)) while !isempty(arr) s = size(arr, 1) - red = dropdims(sum([_dominates(arr[i, :], arr[j, :]) for i in 1:s, j in 1:s], - dims = 1) .== 0, dims = 1) + red = dropdims( + sum([_dominates(arr[i, :], arr[j, :]) for i in 1:s, j in 1:s], + dims = 1) .== 0, + dims = 1) a = 1:s sel::Array{Int64, 1} = a[red] push!(fronts, ind[sel]) @@ -2045,7 +2044,8 @@ function surrogate_optimize(obj, rtea::RTEA, lb, ub, surrRTEAND::AbstractSurroga return pareto_set, pareto_front end -function surrogate_optimize(obj::Function, ::EI, lb::AbstractArray, ub::AbstractArray, krig, +function surrogate_optimize( + obj::Function, ::EI, lb::AbstractArray, ub::AbstractArray, krig, sample_type::SectionSample; maxiters = 100, num_new_samples = 100) dtol = 1e-3 * norm(ub - lb) diff --git a/src/Radials.jl b/src/Radials.jl index 9cff6bf43..f84d67b74 100644 --- a/src/Radials.jl +++ b/src/Radials.jl @@ -42,8 +42,6 @@ of a polynomial term. References: https://en.wikipedia.org/wiki/Polyharmonic_spline - - """ function RadialBasis(x, y, lb, ub; rad::RadialFunction = linearRadial(), scale_factor::Real = 0.5, sparse = false) @@ -112,9 +110,9 @@ using ChainRulesCore: @non_differentiable function _make_combination(n, d, ix) exponents_combinations = [e for e - in collect(Iterators.product(Iterators.repeated(0:n, - d)...))[:] - if sum(e) <= n] + in collect(Iterators.product(Iterators.repeated(0:n, + d)...))[:] + if sum(e) <= n] return exponents_combinations[ix + 1] end @@ -131,12 +129,15 @@ degree `n` and `d` dimensions. Time complexity: `(n+1)^d.` # Example + For n=2, d=2 the multivariate polynomial basis is + ```` 1, x,y x^2,y^2,xy ```` + Therefore the 3rd (ix=3) element is `y` . Therefore when x=(13,43) and ix=3 this function will return 43. """ @@ -145,7 +146,7 @@ function multivar_poly_basis(x, ix, d, n) return one(eltype(x)) else prod(a^d - for (a, d) + for (a, d) in zip(x, _make_combination(n, d, ix))) end end diff --git a/src/Sampling.jl b/src/Sampling.jl index 0eda577e0..6a181b11e 100644 --- a/src/Sampling.jl +++ b/src/Sampling.jl @@ -18,13 +18,14 @@ end #### SectionSample #### """ SectionSample{T}(x0, sa) + `SectionSample(x0, sampler)` where `sampler` is any sampler above and `x0` is a vector of either `NaN` for a free dimension or some scalar for a constrained dimension. """ struct SectionSample{ R <: Real, I <: Integer, VR <: AbstractVector{R}, - VI <: AbstractVector{I}, + VI <: AbstractVector{I} } <: SamplingAlgorithm x0::VR sa::SamplingAlgorithm @@ -36,6 +37,7 @@ free_dimensions(section_sampler::SectionSample)::Vector{Int64} = findall(x -> x isnan.(section_sampler.x0)) """ sample(n,lb,ub,K::SectionSample) + Returns Tuples constrained to a section. In surrogate-based identification and control, optimization can alternate between unconstrained sampling in the full-dimensional parameter space, and sampling constrained on specific sections (e.g. a planes in a 3D volume), A SectionSample allows sampling and optimizing on a subset of 'free' dimensions while keeping 'fixed' ones constrained. @@ -80,6 +82,7 @@ end """ SectionSample(n, d, K::SectionSample) + In surrogate-based identification and control, optimization can alternate between unconstrained sampling in the full-dimensional parameter space, and sampling constrained on specific sections (e.g. planes in a 3D volume). `SectionSample` allows sampling and optimizing on a subset of 'free' dimensions while keeping 'fixed' ones constrained. The sampler is defined diff --git a/src/Surrogates.jl b/src/Surrogates.jl index c9d5a1de1..61e4a62dc 100755 --- a/src/Surrogates.jl +++ b/src/Surrogates.jl @@ -85,22 +85,22 @@ export current_surrogates export GEKPLS export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure export LobachevskyStructure, - NeuralStructure, RandomForestStructure, - SecondOrderPolynomialStructure + NeuralStructure, RandomForestStructure, + SecondOrderPolynomialStructure export WendlandStructure export AbstractSurrogate, SamplingAlgorithm export Kriging, RadialBasis, add_point!, std_error_at_point # Parallelization Strategies export potential_optimal_points export MinimumConstantLiar, MaximumConstantLiar, MeanConstantLiar, KrigingBeliever, - KrigingBelieverUpperBound, KrigingBelieverLowerBound + KrigingBelieverUpperBound, KrigingBelieverLowerBound # radial basis functions export linearRadial, cubicRadial, multiquadricRadial, thinplateRadial # samplers export sample, GridSample, RandomSample, SobolSample, LatinHypercubeSample, - HaltonSample + HaltonSample export RandomSample, KroneckerSample, GoldenSample, SectionSample # Optimization algorithms @@ -112,8 +112,8 @@ export SecondOrderPolynomialSurrogate export Wendland export RadialBasisStructure, KrigingStructure, LinearStructure, InverseDistanceStructure export LobachevskyStructure, - NeuralStructure, RandomForestStructure, - SecondOrderPolynomialStructure + NeuralStructure, RandomForestStructure, + SecondOrderPolynomialStructure export WendlandStructure #export MOE export VariableFidelitySurrogate diff --git a/src/VariableFidelity.jl b/src/VariableFidelity.jl index d90928d38..fc554c35e 100644 --- a/src/VariableFidelity.jl +++ b/src/VariableFidelity.jl @@ -147,8 +147,6 @@ end add_point!(varfid::VariableFidelitySurrogate,x_new,y_new) I expect to add low fidelity data to the surrogate. - - """ function add_point!(varfid::VariableFidelitySurrogate, x_new, y_new) if length(varfid.x[1]) == 1