Add Gaussian Processes

Project.toml

@@ -5,6 +5,7 @@ authors = ["Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>", "Ander Gray <
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918"
 Bootstrap = "e28b5b4c-05e8-5b66-bc03-6f0c0a0a06e0"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 Copulas = "ae264745-0b69-425e-9d9d-cf662c5eec93"
@@ -25,7 +26,9 @@ MeshAdaptiveDirectSearch = "f4d74008-4565-11e9-04bd-4fe404e6a92a"
 Monomials = "272bfe72-f66c-432f-a94d-600f29493792"
 Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
 Mustache = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70"
+NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+ParameterHandling = "2412ca09-6db7-441c-8e3a-88d5709968c5"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
@@ -39,6 +42,7 @@ TransportMaps = "6fd49bf0-03ac-4dd1-b9b3-42e511274d6e"
 
 [compat]
 ADTypes = "1.21.0"
+AbstractGPs = "0.5.24"
 Bootstrap = "2.2"
 Clarabel = "0.11.1"
 Copulas = "0.1.34"
@@ -56,6 +60,7 @@ MeshAdaptiveDirectSearch = "0.1.0"
 Monomials = "1.0"
 Mooncake = "0.4.196, 0.5"
 Mustache = "1.0"
+NLSolversBase = "8.0.1"
 Optim = "2.0, 2.1"
 Primes = "0.5"
 QuadGK = "2.11.1"

demo/metamodels/gaussianprocess.jl

@@ -0,0 +1,42 @@
+using UncertaintyQuantification
+
+x = RandomVariable.(Uniform(-5, 5), [:x1, :x2])
+
+himmelblau = Model(
+    df -> (df.x1 .^ 2 .+ df.x2 .- 11) .^ 2 .+ (df.x1 .+ df.x2 .^ 2 .- 7) .^ 2, :y
+)
+
+design = LatinHypercubeSampling(80)
+
+mean_f = ConstMean(0.0)
+kernel = SqExponentialKernel()
+
+gp_prior = GP(mean_f, kernel)
+
+using Optim
+
+optimizer = MaximumLikelihoodEstimation(Optim.Adam(alpha = 0.005), Optim.Options(; iterations = 10, show_trace = false))
+
+input_transform = ZScoreTransformChoice()
+
+gp_model = GaussianProcess(
+    gp_prior,
+    x,
+    himmelblau,
+    :y;
+    experimental_design = design,
+    input_transform = input_transform,
+    optimizer = optimizer
+)
+
+test_data = sample(x, 1000)
+evaluate!(gp_model, test_data; mode = :mean_and_var)
+
+test_data = sample(x, 1000)
+evaluate!(gp_model, test_data)
+evaluate!(himmelblau, test_data)
+
+mse = mean((test_data.y .- test_data.y_mean) .^ 2)
+println("MSE is:  $mse")
+
+# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

demo/metamodels/gaussianprocessadaptive.jl

@@ -0,0 +1,55 @@
+using UncertaintyQuantification
+using Plots
+using DataFrames
+using Optim # hide
+
+x = RandomVariable.(Uniform(-5, 5), [:x1, :x2])
+himmelblau = Model(
+    df -> (df.x1 .^ 2 .+ df.x2 .- 11) .^ 2 .+ (df.x1 .+ df.x2 .^ 2 .- 7) .^ 2, :y
+)
+
+design = LatinHypercubeSampling(80)
+mean_f = ConstMean(0.0)
+kernel = SqExponentialKernel()
+
+gp_prior = GP(mean_f, kernel)
+input_transform = ZScoreTransformChoice()
+optimizer = MaximumLikelihoodEstimation(Optim.Adam(alpha = 0.005), Optim.Options(; iterations = 10, show_trace = false))
+
+initial_gp = GaussianProcess(
+    gp_prior,
+    x,
+    himmelblau,
+    :y;
+    experimental_design = design,
+    input_transform = input_transform,
+    optimizer = optimizer
+)
+
+learning_function = MaximinDistance()
+n_added_points = 20
+
+adaptive_gp = AdaptiveGaussianProcess(
+    deepcopy(initial_gp),
+    x,
+    himmelblau,
+    learning_function,
+    n_added_points;
+    optimizer = optimizer
+)
+
+test_data = sample(x, LatinHypercubeSampling(1000))
+test_data_adaptive = deepcopy(test_data)
+evaluate!(initial_gp, test_data; mode = :mean)
+evaluate!(himmelblau, test_data)
+
+mse = mean((test_data.y .- test_data.y_mean) .^ 2)
+println("MSE (initial GP):  $mse")
+
+evaluate!(adaptive_gp, test_data_adaptive; mode = :mean)
+evaluate!(himmelblau, test_data_adaptive)
+
+mse_adap = mean((test_data_adaptive.y .- test_data_adaptive.y_mean) .^ 2)
+println("MSE (adaptive GP):  $mse_adap")
+
+# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

docs/Project.toml

@@ -6,6 +6,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 UncertaintyQuantification = "7183a548-a887-11e9-15ce-a56ab60bad7a"
@@ -14,3 +15,4 @@ UncertaintyQuantification = "7183a548-a887-11e9-15ce-a56ab60bad7a"
 Documenter = "1.14.1"
 DocumenterCitations = "1.4.1"
 DocumenterVitepress = "0.2.6"
+Optim = "1.9.4, 2.0"

docs/literate/metamodels/gaussianprocess.jl

@@ -0,0 +1,136 @@
+#===
+# Gaussian Process Regression
+
+## Himmelblau's Function
+
+In this example, we will model the following test function (known as Himmelblau's function) in the range ``x1, x2 ∈ [-5, 5]`` with a Gaussian process (GP) regression model.
+
+It is defined as:
+
+ ```math
+f(x1, x2) = (x1^2 + x2 - 11)^2 + (x1 + x2^2 - 7)^2.
+```
+===#
+# ![](himmelblau.svg)
+#===
+Analogue to the response surface example, we create an array of random variables, that will be used when evaluating the points that our experimental design produces.
+===#
+
+using UncertaintyQuantification
+
+x = RandomVariable.(Uniform(-5, 5), [:x1, :x2])
+
+himmelblau = Model(
+    df -> (df.x1 .^ 2 .+ df.x2 .- 11) .^ 2 .+ (df.x1 .+ df.x2 .^ 2 .- 7) .^ 2, :y
+)
+#md nothing # hide
+
+#===
+Next, we chose a experimental design. In this example, we are using a `LatinHyperCube` design from which we draw 80 samples to train our model:
+===#
+
+design = LatinHypercubeSampling(80)
+
+#===
+After that, we construct a prior GP model. Here we assume a constant mean of 0.0 and a squared exponential kernel with automatic relevance determination (ARD).
+We also assume a small Gaussian noise term in the observations for numerical stability:
+===#
+
+mean_f = ConstMean(0.0)
+kernel = SqExponentialKernel()
+
+gp_prior = GP(mean_f, kernel)
+
+#===
+Next, we set up an optimizer used in the log marginal likelihood maximization to find the optimal hyperparameters of our GP model. Here we use the Adam optimizer from the `Optim.jl` package with a learning rate of 0.005 and run it for 10 iterations.:
+===#
+using Optim
+
+optimizer = MaximumLikelihoodEstimation(Optim.Adam(alpha = 0.005), Optim.Options(; iterations = 10, show_trace = false))
+#md nothing # hide
+
+#===
+Finally, we define an input standardization (here a z-score transform). While not strictly necessary for this example, standardization can help finding good hyperparameters.
+Note that we can also define an output transform to scale the output for training the GP. When evaluating the GP model, the input will be automatically transformed with the fitted standardization.
+The output will be transformed back to the original scale automatically as well.
+===#
+
+input_transform = ZScoreTransformChoice()
+#md nothing # hide
+
+#===
+The GP regression model is now constructed by calling the `GaussianProcess` constructor with the prior GP, the input random variables, the model, the output symbol, the experimental design, and the optional input and output transform choices.
+The construction then samples the experimental design, evaluates the model at the sampled points, standardizes the input and output data, and constructs the posterior GP.
+===#
+#md using Random #hide
+#md Random.seed!(42) #hide
+
+gp_model = GaussianProcess(
+    gp_prior,
+    x,
+    himmelblau,
+    :y;
+    experimental_design = design,
+    input_transform = input_transform,
+    optimizer = optimizer
+)
+#md nothing # hide
+
+#===
+The GP regression model uses finite projections of the fitted posterior GP to make predictions. As of now, the hyperparameters of the GP might not be optimal.
+We can find optimal hyperparameters through maximizing the log marginal likelihood of observing the training data under the posterior GP.
+===#
+
+#===
+To evaluate the `GaussianProcess`, use `evaluate!(gp::GaussianProcess, data::DataFrame)` with the `DataFrame` containing the points you want to evaluate.
+The evaluation of a GP is not unique, and we can choose to evaluate the mean prediction, the prediction variance, a combination of both, or draw samples from the posterior distribution.
+The default is to evaluate the mean prediction.
+We can specify the evaluation mode via the `mode` keyword argument. Supported options are:
+- `:mean` - predictive mean (default)
+- `:var` - predictive variance
+- `:mean_and_var` - both mean and variance
+- `:sample` - random samples from the predictive distribution
+===#
+
+test_data = sample(x, 1000)
+evaluate!(gp_model, test_data; mode = :mean_and_var)
+
+#===
+The mean prediction of our model in this case has an mse of about 65 and looks like this in comparison to the original:
+===#
+
+#md using Plots #hide
+#md using DataFrames #hide
+#md a = range(-5, 5; length=200) #hide
+#md b = range(-5, 5; length=200) #hide
+#md A = repeat(collect(a)', length(b), 1) #hide
+#md B = repeat(collect(b), 1, length(a)) #hide
+#md df = DataFrame(x1 = vec(A), x2 = vec(B)) #hide
+#md evaluate!(gp_model, df; mode=:mean_and_var) #hide
+#md evaluate!(himmelblau, df) #hide
+#md gp_mean = reshape(df[:, :y_mean], length(b), length(a)) #hide
+#md gp_var = reshape(df[:, :y_var], length(b), length(a)) #hide
+#md himmelblau_values = reshape(df[:, :y], length(b), length(a)) #hide
+#md s1 = surface(a, b, himmelblau_values; plot_title="Himmelblau's function")
+#md s2 = surface(a, b, gp_mean; plot_title="GP posterior mean")
+#md plot(s1, s2, layout = (1, 2), legend = false)
+#md savefig("gp-mean-comparison.svg") # hide
+#md s3 = surface(a, b, gp_var; plot_title="GP posterior variance") # hide
+#md plot(s3, legend = false) #hide
+#md savefig("gp-variance.svg"); nothing # hide
+
+# ![](gp-mean-comparison.svg)
+
+#===
+Note that the mse in comparison to the response surface model (with an mse of about 1e-26) is significantly higher.
+However, the GP model also provides a measure of uncertainty in its predictions via the predictive variance.
+===#
+
+# ![](gp-variance.svg)
+
+#jl test_data = sample(x, 1000)
+#jl evaluate!(gp_model, test_data)
+#jl evaluate!(himmelblau, test_data)
+
+#jl mse = mean((test_data.y .- test_data.y_mean) .^ 2)
+#jl println("MSE is:  $mse")

docs/literate/metamodels/gaussianprocessadaptive.jl

@@ -0,0 +1,105 @@
+#===
+# Adaptive Gaussian Process Regression
+
+Adaptive Gaussian process regression enriches an initial surrogate model with
+new evaluations selected by a learning function. At every iteration, the
+algorithm samples candidate points, selects the most informative one according
+to the learning function, evaluates the expensive model there, and refits the
+Gaussian process.
+
+## Himmelblau's Function
+
+As in the (non-adaptive) GP example, we consider the Himmelblau function in ``x1, x2 ∈ [-5, 5]``
+as a test function.
+===#
+
+#md using UncertaintyQuantification # hide
+#md using Plots # hide
+#md using DataFrames # hide
+#md using Optim # hide
+
+#jl using UncertaintyQuantification
+#jl using Plots
+#jl using DataFrames
+#jl using Optim # hide
+
+#===
+First, define the probabilistic input and the expensive model to approximate.
+===#
+
+x = RandomVariable.(Uniform(-5, 5), [:x1, :x2])
+himmelblau = Model(
+    df -> (df.x1 .^ 2 .+ df.x2 .- 11) .^ 2 .+ (df.x1 .+ df.x2 .^ 2 .- 7) .^ 2, :y
+)
+#md nothing # hide
+
+#===
+We start with the same initial Gaussian process surrogate as in the *regular* GP regression
+example. Hence, we use the same initial design and same optimizer.
+===#
+
+design = LatinHypercubeSampling(80)
+mean_f = ConstMean(0.0)
+kernel = SqExponentialKernel()
+
+gp_prior = GP(mean_f, kernel)
+input_transform = ZScoreTransformChoice()
+optimizer = MaximumLikelihoodEstimation(Optim.Adam(alpha = 0.005), Optim.Options(; iterations = 10, show_trace = false))
+
+initial_gp = GaussianProcess(
+    gp_prior,
+    x,
+    himmelblau,
+    :y;
+    experimental_design = design,
+    input_transform = input_transform,
+    optimizer = optimizer
+)
+#md nothing # hide
+
+#===
+Next, we update the initial GP using a selected learning function and a set number of
+additional points used to refine the initial GP. We use the [`MaximinDistance`](@ref)
+acquisition function and select `20` additional points.
+===#
+
+learning_function = MaximinDistance()
+n_added_points = 20
+#md nothing # hide
+
+#===
+We refine the GP using [`AdaptiveGaussianProcess`](@ref) which we pass our initial GP, the
+`learning_function` and `n_added_points`.
+===#
+
+adaptive_gp = AdaptiveGaussianProcess(
+    deepcopy(initial_gp),
+    x,
+    himmelblau,
+    learning_function,
+    n_added_points;
+    optimizer = optimizer
+)
+#md nothing # hide
+
+#===
+To assess the fitted surrogate, we compute the MSE between GP mean and the reference model.
+We compare the MSE of the initial GP and the refined GP.
+
+We start with the initial GP:
+===#
+
+test_data = sample(x, LatinHypercubeSampling(1000))
+test_data_adaptive = deepcopy(test_data)
+evaluate!(initial_gp, test_data; mode = :mean)
+evaluate!(himmelblau, test_data)
+
+mse = mean((test_data.y .- test_data.y_mean) .^ 2)
+println("MSE (initial GP):  $mse")
+
+# Then, we also evaluate the adaptively refined GP at the same test set:
+evaluate!(adaptive_gp, test_data_adaptive; mode = :mean)
+evaluate!(himmelblau, test_data_adaptive)
+
+mse_adap = mean((test_data_adaptive.y .- test_data_adaptive.y_mean) .^ 2)
+println("MSE (adaptive GP):  $mse_adap")

docs/make.jl

@@ -58,6 +58,7 @@ makedocs(;
             "Basis Functions" => "api/basisfunctions.md",
             "ResponseSurface" => "api/responsesurface.md",
             "PolyharmonicSpline" => "api/polyharmonicspline.md",
+            "Gaussian Processes" => "api/gaussianprocesses.md",
             "Simulations" => "api/simulations.md",
             "Bayesian Updating" => "api/bayesianupdating.md",
             "Transport Maps" => "api/transportmaps.md",