diff --git a/Project.toml b/Project.toml index 148999127..59d060175 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ authors = ["Jasper Behrensdorf ", "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" diff --git a/demo/metamodels/gaussianprocess.jl b/demo/metamodels/gaussianprocess.jl new file mode 100644 index 000000000..b919ada46 --- /dev/null +++ b/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 diff --git a/demo/metamodels/gaussianprocessadaptive.jl b/demo/metamodels/gaussianprocessadaptive.jl new file mode 100644 index 000000000..db1d91cfd --- /dev/null +++ b/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 diff --git a/docs/Project.toml b/docs/Project.toml index ac218a511..165e8ba07 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/literate/metamodels/gaussianprocess.jl b/docs/literate/metamodels/gaussianprocess.jl new file mode 100644 index 000000000..90bbcfb88 --- /dev/null +++ b/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") diff --git a/docs/literate/metamodels/gaussianprocessadaptive.jl b/docs/literate/metamodels/gaussianprocessadaptive.jl new file mode 100644 index 000000000..48147af72 --- /dev/null +++ b/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") diff --git a/docs/make.jl b/docs/make.jl index 3811c96b4..7891de6e5 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/references.bib b/docs/references.bib index 005b59131..d686722e4 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -74,6 +74,20 @@ @article{biBhattacharyyaDistanceEnriching2019 doi = {10.1016/j.ymssp.2019.04.035} } +@article{bichon2008efficient, + title = {Efficient {{Global Reliability Analysis}} for {{Nonlinear Implicit Performance Functions}}}, + author = {Bichon, B. J. and Eldred, M. S. and Swiler, L. P. and Mahadevan, S. and McFarland, J. M.}, + year = {2008}, + month = oct, + journal = {AIAA Journal}, + volume = {46}, + number = {10}, + pages = {2459--2468}, + issn = {0001-1452, 1533-385X}, + doi = {10.2514/1.34321}, +} + + @article{burkEfficientSampling, author = {Burk, Kyle M. and Narayan, Akil and Orr, Joseph A.}, title = {Efficient sampling for polynomial chaos-based uncertainty quantification and sensitivity analysis using weighted approximate Fekete points}, @@ -96,9 +110,6 @@ @article{chan2022adaptive publisher = {Elsevier} } - - - @article{chingTransitionalMarkovChain2007, title = {Transitional {{Markov Chain Monte Carlo Method}} for {{Bayesian Model Updating}}, {{Model Class Selection}}, and {{Model Averaging}}}, author = {Ching, Jianye and Chen, Yi-Chu}, @@ -169,6 +180,20 @@ @article{dubois1990consonant publisher = {Elsevier} } +@article{echard2011akmcs, + title = {{{AK-MCS}}: {{An}} Active Learning Reliability Method Combining {{Kriging}} and {{Monte Carlo Simulation}}}, + author = {Echard, B. and Gayton, N. and Lemaire, M.}, + year = {2011}, + month = mar, + journal = {Structural Safety}, + volume = {33}, + number = {2}, + pages = {145--154}, + issn = {01674730}, + doi = {10.1016/j.strusafe.2011.01.002}, + urldate = {2023-10-25}, +} + @techreport{fersonConstructingProbabilityBoxes2015, title = {Constructing {{Probability Boxes}} and {{Dempster-Shafer Structures}}}, author = {Ferson, Scott and Kreinovick, Vladik and Ginzburg, Lev and Sentz, Fari and Meyers, Davis}, @@ -179,6 +204,20 @@ @techreport{fersonConstructingProbabilityBoxes2015 doi = {10.2172/809606} } +@article{fuhg2021state, + title = {State-of-the-{{Art}} and {{Comparative Review}} of {{Adaptive Sampling Methods}} for {{Kriging}}}, + author = {Fuhg, Jan N. and Fau, Am{\'e}lie and Nackenhorst, Udo}, + year = 2021, + month = jun, + journal = {Archives of Computational Methods in Engineering}, + volume = {28}, + number = {4}, + pages = {2689--2747}, + issn = {1886-1784}, + doi = {10.1007/s11831-020-09474-6}, +} + + @phdthesis{grashornEfficientDiagnostics2024, type = {{{Doctoral Thesis}}}, title = {Efficient Diagnostics of Complex Mechanical Systems}, @@ -454,6 +493,20 @@ @book{raiffaAppliedStatisticalDecision1961 pagetotal = {356} } +@book{rasmussen2005gaussian, + title = {Gaussian {Processes} for {Machine} {Learning}}, + copyright = {http://creativecommons.org/licenses/by-nc-nd/4.0/}, + isbn = {978-0-262-25683-4}, + url = {https://direct.mit.edu/books/book/2320/Gaussian-Processes-for-Machine-Learning}, + language = {en}, + urldate = {2025-09-04}, + publisher = {The MIT Press}, + author = {Rasmussen, Carl Edward and Williams, Christopher K. I.}, + month = nov, + year = {2005}, + doi = {10.7551/mitpress/3206.001.0001} +} + @misc{ramgraberTriangularTransport2025, title = {A Friendly Introduction to Triangular Transport}, author = {Ramgraber, Maximilian and Sharp, Daniel and Provost, Mathieu Le and Marzouk, Youssef}, diff --git a/docs/src/api/gaussianprocesses.md b/docs/src/api/gaussianprocesses.md new file mode 100644 index 000000000..62e060319 --- /dev/null +++ b/docs/src/api/gaussianprocesses.md @@ -0,0 +1,35 @@ +# Gaussian Process Regression + +Methods for Gaussian process regression. + +## Index + +```@index +Pages = ["gaussianprocesses.md"] +``` + +## Types + +```@docs +GaussianProcess +MaximumLikelihoodEstimation +IdentityTransformChoice +ZScoreTransformChoice +UnitRangeTransformChoice +StandardNormalTransformChoice +MaximumVariance +ExpectedImprovement +ProbabilityOfImprovement +UpperConfidenceBound +DeviationNumber +ExpectedFeasibility +MaximinDistance +ExpectedImprovementForGlobalFit +``` + +## Functions + +```@docs +AdaptiveGaussianProcess +evaluate!(gp::GaussianProcess, data::DataFrame; mode::Symbol = :mean, n_samples::Int = 1) +``` diff --git a/docs/src/manual/metamodels.md b/docs/src/manual/metamodels.md index b1a153e8b..983dfc822 100644 --- a/docs/src/manual/metamodels.md +++ b/docs/src/manual/metamodels.md @@ -235,3 +235,366 @@ The reliability of the IPM, that is the probability that and unobserved data poi ### Reliability analysis As the IPM is an imprecise model, it can only be applied in a reliability analysis using the [`DoubleLoop`](@ref) or [`RandomSlicing`](@ref). For more information, see [Imprecise Reliability Analysis](@ref). + +## Gaussian Process Regression + +### Theoretical Background + +A Gaussian Process (GP) is a collection of random variables, any finite subset of which has a joint Gaussian distribution. It is fully specified by a mean function $m(x)$ and a covariance (kernel) function $k(x, x')$. In GP regression, we aim to model an unknown function $f(x)$. Before observing any data, we assume that the function $f(x)$ is distributed according to a GP: + +```math +f(x) \sim \mathcal{G}\mathcal{P}\left( m(x), k(x, x') \right). +``` + +This prior GP specifies that any finite collection of function values follows a multivariate normal distribution. + +To define a prior GP we use [`AbstractGPs.jl`](https://juliagaussianprocesses.github.io/AbstractGPs.jl/stable/) for the GP interface and mean function, and [`KernelFunctions.jl`](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/) for the definition of a covariance kernel. Below, we construct a simple prior GP with a constant zero mean function and a scaled squared exponential kernel: + +```@example gaussianprocess +using UncertaintyQuantification + +kernel = SqExponentialKernel() ∘ ScaleTransform(3.0) +gp = GP(0.0, kernel); nothing # hide +``` + +Note that the definition of a prior GP is handled by `UncertaintyQuantification` if no prior GP is specified. The construction of a `GaussianProcess` is flexible. Mean functions, kernels and many other parameters can be specified later directly in the constructor of the `GaussianProcess`. + +#### Posterior Gaussian Process + +The posterior GP represents the distribution of functions after incorporating observed data. We denote the observation data as: + +```math +\mathcal{D} = \lbrace (\hat{x}_i, \hat{f}_i) \mid i=1, \dots, N \rbrace, +``` + +where $\hat{f}_i = f(\hat{x}_i)$ in the noise-free observation case, and $\hat{f}_i = f(\hat{x}_i) + e_i$ in the noisy case, with independent noise terms $e_i \sim \mathcal{N}(0, \sigma_e^2)$. Let $\hat{X} = [\hat{x}_1, \dots, \hat{x}_N]$ denote the collection of observation data locations. The corresponding mean vector and covariance matrix are: + +```math +\mu(\hat{X}) = [m(\hat{x}_1), \dots, m(\hat{x}_N)], \quad K(\hat{X}, \hat{X}) \text{ with entries } K_{ij} = k(\hat{x}_i, \hat{x}_j). +``` + +For a new input location $x^*$ we are interested at the unknown function value $f^* = f(x^*)$. By the definition of a GP, the joint distribution of observed outputs $\hat{f}_i$ and the unknown $f^*$ is multivariate Gaussian: + +```math +\begin{bmatrix} \hat{f}\\ f^* \end{bmatrix} = \mathcal{N}\left( \begin{bmatrix} \mu(\hat{X}) \\ m(x^*) \end{bmatrix}, \begin{bmatrix} K(\hat{X}, \hat{X}) & K(\hat{X}, x^*)\\ K(x^*, \hat{X}) & K(x^*, x^*) \end{bmatrix} \right), +``` + +where: + +- ``K(\hat{X}, \hat{X})`` is the covariance matrix with entries ``K_{ij} = k(\hat{x}_i, \hat{x}_j)``, +- ``K(\hat{X}, x^*)`` is the covariance matrix with entries ``K_{i1} = k(\hat{x}_i, x^*)``, +- and ``(x^*, x^*)`` is the variance at the unknown input location. + +We can then obtain the posterior distribution of $f^*$ from the properties of multivariate Gaussian distributions (see, e.g. Appendix A.2 in [rasmussen2005gaussian](@cite)), by conditioning the joint Gaussian on the observed outputs $\hat{f}_i$: + +```math +f^* \mid \hat{X}, \hat{f}, x^* \sim \mathcal{N}(\mu^*(x^*), \Sigma^*(x^*)), +``` + +with + +```math +\mu^*(x^*) = m(x^*) + K(x^*, \hat{X})K(\hat{X}, \hat{X})^{-1}(\hat{f} - \mu(\hat{X})), \\ +\Sigma^*(x^*) = K(x^*, x^*) - K(x^*, \hat{X})K(\hat{X}, \hat{X})^{-1}K(\hat{X}, x^*). +``` + +In the noisy observation case, the covariance between training points is adjusted by adding the noise variance: + +```math +K(\hat{X}, \hat{X}) \rightarrow K(\hat{X}, \hat{X}) + \sigma^2_{e}I. +``` + +The computation of the posterior predictive distribution generalizes straightforwardly to multiple input locations, providing both the posterior mean, which can serve as a regression estimate of the unknown function, and the posterior variances, which quantify the uncertainty at each point. Because the posterior is multivariate Gaussian, one can also sample function realizations at specified locations to visualize possible functions consistent with the observed data. + +To construct a posterior GP, we need to define training data in form of a `DataFrame`. Constructing a `GaussianProcess` model will then automatically compute the posterior GP to predict requested the modeled output $y$ and by default it will also optimize the hyperparameters. If this is not desired, the input `learn_hyperparameters=false` can be set. + +The following creates a standard GP with mean function `ConstMean()`, kernel `SqExponentialKernel()`, and directly optimizes the hyperparameters. Note that while `ConstMean(0.0)` and `ZeroMean()` provide the same zero-mean prior GP, using `ConstMean()` also allows for optimization of the mean. +We also equip the GP with small observation noise $\sigma^2$, which has implications on the numerical stability and allows the GP to handle imprecise data. The noise can also be optimized as part of the hyperparameter optimization, but it is not optimized by default. +To specify different mean functions and/or kernels, either construct a GP manually beforehand, or use them as inputs. + +```@example gaussianprocess +using DataFrames # hide +x = collect(range(0, 10, 10)) +y = sin.(x) + 0.3 * cos.(2 .* x) +df = DataFrame(x = x, y = y) + +mean_fct = ConstMean(0.0) +kernel = SqExponentialKernel() ∘ ScaleTransform(3.0) + +gp_prior = GP(mean_fct, kernel) + +σ² = 1e-5 + +# these are equivalent +gp_model = GaussianProcess(gp_prior, df, :y; σ²=σ²) +gp_model = GaussianProcess(df, :y; σ²=σ², mean_fct=mean_fct, kernel=kernel) +# providing the input learn_noise=true also optimizes the data noise +gp_model = GaussianProcess(df, :y; σ²=σ², mean_fct=mean_fct, kernel=kernel, learn_noise=true); nothing # hide +``` + +Now we can use our GP model to predict at new input locations `x_test`: + +```@example gaussianprocess +using Plots # hide +x_test = collect(range(0, 5, 500)) +prediction = DataFrame(:x => x_test) + +evaluate!(gp_model, prediction; mode=:mean_and_var) + +prediction_mean = prediction[!, :y_mean] # hide +prediction_std = sqrt.(prediction[!, :y_var]) # hide + +p = plot(x_test, prediction_mean, color=:blue, label="Mean prediction") # hide +plot!( + x_test, prediction_mean, ribbon=2 .* prediction_std, + color=:grey, alpha=0.5, label="Confidence band" +) # hide + +y_true = sin.(x_test) + 0.3 * cos.(2 .* x_test) # hide +plot!(x_test, y_true, color=:red, label="True function") # hide + +savefig(p, "posterior-gp.svg"); nothing # hide +``` + +![Fitted Gaussian process](posterior-gp.svg) + +#### Hyperparameter optimization + +GP models typically contain hyperparameters in their mean functions $m(x; \theta_m)$ and covariance kernel functions $k(x, x'; \theta_k)$. The observation noise variance $\sigma^2_{e}$ is also considered a hyperparameter related to the kernel. The choice of hyperparameters strongly affects the quality of the posterior GP. + +A common approach to selecting hyperparameters is maximum likelihood estimation (MLE) (see, e.g. [rasmussen2005gaussian](@cite)), where we maximize the likelihood of observing the training data $\mathcal{D}$ under the chosen GP prior. + +The marginal likelihood of the observed training outputs $\hat{f}$ is: + +```math +p(\hat{f} \mid \hat{X}, \theta_m, \theta_k, \sigma^2_{e}) = \mathcal{N}(\hat{f} \mid \mu_{\theta_m}(\hat{X}), K_{\theta_k}(\hat{X}, \hat{X}) + \sigma^2_{e}I), +``` + +where $\mu_{\theta_m}(\hat{X})$ and $K_{\theta_k}(\hat{X}, \hat{X})$ denote the parameter dependent versions of the previously defined quantities. + +For numerical reasons, the logarithm of the marginal likelihood is typically used. Maximizing the log marginal likelihood with respect to the hyperparameters then yields the parameters that best explain the observed data. After obtaining the optimal hyperparameters, the posterior GP can be constructed as described above. + +`UncertaintyQuantification.jl` provides a default optimizer for the hyperparameters based on the [`MaximumLikelihoodEstimation`](@ref) constructor. + +``` +optimizer::AbstractHyperparameterOptimization=MaximumLikelihoodEstimation(Optim.LBFGS(), Optim.Options(; iterations=100, show_trace=false)) +``` + +If other options are desired, a different optimizer can be constructed based on [`Optim.jl`](https://julianlsolvers.github.io/Optim.jl/stable/). The script below shows the difference between an optimized and unoptimized GP. + +```@example gaussianprocess +using Optim + +optimization = MaximumLikelihoodEstimation( + Optim.LBFGS(), + Optim.Options(; iterations=10, show_trace=false) + ) + +gp_model = GaussianProcess(df, :y; + σ²=σ², + mean_fct=mean_fct, + kernel=kernel, + optimizer=optimization + ) + +gp_model_unoptimized = GaussianProcess(df, :y; + σ²=σ², + mean_fct=mean_fct, + kernel=kernel, + learn_hyperparameters=false + ) + +prediction = DataFrame(:x => x_test) +prediction_unopt = DataFrame(:x => x_test) +evaluate!(gp_model, prediction; mode=:mean_and_var) +evaluate!(gp_model_unoptimized, prediction_unopt; mode=:mean_and_var) + +prediction_mean = prediction[!, :y_mean] # hide +prediction_std = sqrt.(prediction[!, :y_var]) # hide + +prediction_mean_unopt = prediction_unopt[!, :y_mean] #hide +prediction_std_unopt = sqrt.(prediction_unopt[!, :y_var]) # hide + +p = plot(x_test, prediction_mean, ribbon=2 .* prediction_std, color=:blue, alpha=0.5, label="Optimized") # hide +plot!(x_test, prediction_mean_unopt, ribbon=2 .* prediction_std_unopt, color=:grey, alpha=0.5, label="Not Optimized") # hide + +plot!(x_test, y_true, color=:red, label="True function") # hide + +savefig(p, "posterior-gp-opt.svg"); nothing # hide +``` + +![Optimized Gaussian process](posterior-gp-opt.svg) + +Internally, `MaximumLikelihoodEstimation()` defaults to using [`LBFGS`](https://julianlsolvers.github.io/Optim.jl/stable/algo/lbfgs/) optimizer that performs 100 optimization steps with standard optimization hyperparameters as defined [`Optim.jl`](https://julianlsolvers.github.io/Optim.jl/stable/). Note that any other first-order optimizer supported by [`Optim.jl`](https://julianlsolvers.github.io/Optim.jl/stable/), along with its corresponding hyperparameters, can also be used when constructing [`MaximumLikelihoodEstimation`](@ref). + +During optimization, GP hyperparameters $\theta_m, \theta_k$ and $\sigma^2_{e}$ are automatically extracted and updated. + +We support the automatic extraction of hyperparameters from mean functions provided by [`AbstractGPs.jl`](https://juliagaussianprocesses.github.io/AbstractGPs.jl/stable/api/#Mean-functions), with the exception of: + +- Custom mean functions [`CustomMean`](https://juliagaussianprocesses.github.io/AbstractGPs.jl/stable/api/#AbstractGPs.CustomMean). These are defined with a custom function that itself could depend on hyperparameters. These additional hyperparameters are ignored in the optimization. + +Kernel functions are defined with the kernels and transformations provided by [`KernelFunctions.jl`](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/). For similar reasons as with `CustomMean`, we do not extract potential function hyperparameters from the following kernels or transforms: + +- Transforms defined with custom functions [`FunctionTransform`](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/transform/#KernelFunctions.FunctionTransform), +- The [`GibbsKernel`](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/kernels/#KernelFunctions.GibbsKernel), which models a kernel lengthscale parameter with the help of a function. + +Further, GP models containing the following kernels are not supported for hyperparameter optimization currently: + +- Multi-output kernels [`MOKernel`](https://juliagaussianprocesses.github.io/KernelFunctions.jl/stable/kernels/#Multi-output-Kernels), +- Neural kernel networks [`NeuralKernelNetwork`]. + +## Adaptive Gaussian Process Regression + +Fitting a good GP surrogate can require many expensive model evaluations if the initial +experimental design is chosen naively. **Adaptive** (or *active learning*) Gaussian process +regression instead starts from a small initial design and iteratively enriches the training +data: at each iteration a set of candidate points is sampled from the input space, an +**acquisition function** (also called a *learning function*) scores every candidate, the most +promising candidate is evaluated with the true (expensive) model, and the GP is refitted with +the enlarged training set. This is repeated for a fixed number of iterations, or until the +acquisition function's own convergence criterion is met. + +The [`AdaptiveGaussianProcess`](@ref) function drives this loop. It first constructs (or +accepts) an initial [`GaussianProcess`](@ref), then calls `evaluate!` on the supplied `model` +for each newly selected point. + +```@example adaptivegp +using UncertaintyQuantification # hide + +x = RandomVariable(Uniform(-10, 10), :x1) +model = Model(df -> sin.(df.x1) .* df.x1 .^ 2, :y) + +mean_f = ConstMean(0.0) +kernel = Matern52Kernel() +gp_prior = GP(mean_f, kernel) + +n_design_points = 10 +n_added_points = 5 + +adaptive_gp = AdaptiveGaussianProcess( + gp_prior, + x, + model, + :y, + MaximumVariance(), + n_added_points, + n_design_points, +) +nothing # hide +``` + +As with [`GaussianProcess`](@ref), the initial `n_design_points` are sampled with an +`experimental_design` (`LatinHypercubeSampling` by default), while the `n_added_points` +adaptively selected candidates are drawn from `candidate_sampling`, a Monte Carlo sampling +scheme (`MonteCarlo(100_000)` by default). Hyperparameters can be re-optimized after every +added point via `learn_hyperparameters` (default `true`). + +The resulting `adaptive_gp` is a regular [`GaussianProcess`](@ref) and can be evaluated as usual: + +```@example adaptivegp +using DataFrames +using Plots + +test_data = DataFrame(x1 = -10:0.1:10) +evaluate!(adaptive_gp, test_data; mode = :mean_and_var) +evaluate!(model, test_data) + +p = plot(test_data.x1, test_data.y_mean; ribbon = 2 .* sqrt.(test_data.y_var), label = "GP mean ± 2σ", xlabel = "x₁", ylabel = "y", color = :blue, alpha = 0.5) +plot!(p, test_data.x1, test_data.y; label = "True function", color = :red, linestyle = :dash) +scatter!(p, adaptive_gp.training_data.x1[1:n_design_points], adaptive_gp.training_data.y[1:n_design_points]; label = "Initial design", color = :black) +scatter!(p, adaptive_gp.training_data.x1[(n_design_points + 1):end], adaptive_gp.training_data.y[(n_design_points + 1):end]; label = "Adaptively added") +savefig(p, "adaptive_gp_example.svg"); nothing # hide +``` + +# ![Adaptive GP](adaptive_gp_example.svg) + +### Acquisition Functions + +The acquisition function determines *which* candidate point is added next, and therefore what +the adaptive scheme optimizes for. Given the posterior mean ``\mu(\mathbf{x})`` and posterior +standard deviation ``\sigma(\mathbf{x})`` of the current GP, the next point ``\mathbf{x}^{+}`` +is chosen from the sampled candidates ``\mathbf{x} \in \mathcal{X}_c`` by maximizing (or +minimizing) a criterion ``a(\mathbf{x})``. For a review of various acquisition functions we +refer to [fuhg2021state](@cite). + +`UncertaintyQuantification.jl` provides the following acquisition functions to adapatively +refine GP regression models; we separate them by there main area of application: + +### General active learning +Goal: improve the global fit of the GP + +- [`MaximumVariance`](@ref) simply adds the point of maximum posterior variance, + +```math +\mathbf{x}^{+} = \underset{\mathbf{x} \in \mathcal{X}_c}{\operatorname{argmax}}\ \sigma^2(\mathbf{x}). +``` + +- [`MaximinDistance`](@ref) is a space-filling criterion that adds the candidate farthest (in input space) from every existing training point ``\hat{\mathbf{x}}_i \in \hat{X}``, + +```math +\mathbf{x}^{+} = \underset{\mathbf{x} \in \mathcal{X}_c}{\operatorname{argmax}}\ \min_{i} \lVert \mathbf{x} - \hat{\mathbf{x}}_i \rVert. +``` + +- [`ExpectedImprovementForGlobalFit`](@ref) (EIGF) trades off the local discrepancy to the nearest training observation ``\hat{f}_{i(\mathbf{x})}`` (with ``i(\mathbf{x}) = \operatorname{argmin}_i \lVert \mathbf{x} - \hat{\mathbf{x}}_i \rVert``) against the posterior variance, + +```math +\mathbf{x}^{+} = \underset{\mathbf{x} \in \mathcal{X}_c}{\operatorname{argmax}}\ \left[\mu(\mathbf{x}) - \hat{f}_{i(\mathbf{x})}\right]^2 + \sigma^2(\mathbf{x}). +``` + +### Bayesian optimization +Goal: refine the global minimum ``\hat{f}_{\text{best}} = \min_i \hat{f}_i`` + +- [`ExpectedImprovement`](@ref) with exploration parameter ``\xi``, + +```math +\mathrm{EI}(\mathbf{x}) = \left(\hat{f}_{\text{best}} - \mu(\mathbf{x}) - \xi\right) \Phi(z) + \sigma(\mathbf{x}) \phi(z), \qquad +z = \frac{\hat{f}_{\text{best}} - \mu(\mathbf{x}) - \xi}{\sigma(\mathbf{x})}, +``` + +```math +\mathbf{x}^{+} = \underset{\mathbf{x} \in \mathcal{X}_c}{\operatorname{argmax}}\ \mathrm{EI}(\mathbf{x}), +``` + +where ``\Phi`` and ``\phi`` are the standard normal cdf and pdf, respectively (``\mathrm{EI}(\mathbf{x}) = 0`` if ``\sigma(\mathbf{x}) = 0``). + +- [`ProbabilityOfImprovement`](@ref), + +```math +\mathrm{PI}(\mathbf{x}) = \Phi(z), \qquad z = \frac{\hat{f}_{\text{best}} - \mu(\mathbf{x}) - \xi}{\sigma(\mathbf{x})}, \qquad +\mathbf{x}^{+} = \underset{\mathbf{x} \in \mathcal{X}_c}{\operatorname{argmax}}\ \mathrm{PI}(\mathbf{x}). +``` + +- [`UpperConfidenceBound`](@ref) with exploration weight ``\kappa`` minimizes a lower confidence bound (for a minimization objective), + +```math +\mathbf{x}^{+} = \underset{\mathbf{x} \in \mathcal{X}_c}{\operatorname{argmin}}\ \mu(\mathbf{x}) - \kappa\, \sigma(\mathbf{x}). +``` + +### Reliability analysis +Goal: refine the limit-state surface ``g(\mathbf{x}) = \tau`` (typically, ``\tau = 0``): + +- [`DeviationNumber`](@ref), the ``U``-function used in AK-MCS, adds the point closest to the limit state relative to its uncertainty, + +```math +\mathbf{x}^{+} = \underset{\mathbf{x} \in \mathcal{X}_c}{\operatorname{argmin}}\ U(\mathbf{x}), \qquad +U(\mathbf{x}) = \frac{|\mu(\mathbf{x}) - \tau|}{\sigma(\mathbf{x})}. +``` + +- [`ExpectedFeasibility`](@ref) (EFF) integrates the probability that the true response lies within an ``\epsilon``-band ``\epsilon(\mathbf{x}) = \text{epsilon\_factor} \cdot \sigma(\mathbf{x})`` around the limit state, + +```math +\begin{aligned} +\mathrm{EFF}(\mathbf{x}) = {} & \left(\mu(\mathbf{x}) - \tau\right) \Big[2\Phi(z) - \Phi(z^-) - \Phi(z^+)\Big] \\ +& - \sigma(\mathbf{x}) \Big[2\phi(z) - \phi(z^-) - \phi(z^+)\Big] + \epsilon(\mathbf{x}) \Big[\Phi(z^+) - \Phi(z^-)\Big], +\end{aligned} +``` + +```math +z = \frac{\tau - \mu(\mathbf{x})}{\sigma(\mathbf{x})}, \quad +z^- = \frac{\tau - \epsilon(\mathbf{x}) - \mu(\mathbf{x})}{\sigma(\mathbf{x})}, \quad +z^+ = \frac{\tau + \epsilon(\mathbf{x}) - \mu(\mathbf{x})}{\sigma(\mathbf{x})}, +``` + +```math +\mathbf{x}^{+} = \underset{\mathbf{x} \in \mathcal{X}_c}{\operatorname{argmax}}\ \mathrm{EFF}(\mathbf{x}). +``` diff --git a/src/UncertaintyQuantification.jl b/src/UncertaintyQuantification.jl index 17b9a3eea..6ade647c2 100644 --- a/src/UncertaintyQuantification.jl +++ b/src/UncertaintyQuantification.jl @@ -21,6 +21,7 @@ using Monomials using Mooncake: Mooncake using Mustache using Optim +using ParameterHandling using Primes using QuadGK using QuasiMonteCarlo @@ -31,6 +32,9 @@ using StatsBase using TransportMaps using RecipesBase +import NLSolversBase + +@reexport using AbstractGPs @reexport using TransportMaps @reexport using Distributions @reexport using DifferentiationInterface @@ -85,7 +89,12 @@ abstract type AbstractHPCScheduler end abstract type AbstractTransportMap <: ContinuousMultivariateDistribution end +abstract type AbstractAcquisitionFunction end +abstract type AbstractGaussianProcessAcquisitionFunction <: AbstractAcquisitionFunction end + # Types +export AbstractAcquisitionFunction +export AbstractGaussianProcessAcquisitionFunction export AbstractBayesianMethod export AbstractBayesianPointEstimate export AbstractDesignOfExperiments @@ -103,6 +112,7 @@ export UQModel # Structs export AdvancedLineSampling +export AdaptiveGaussianProcess export EmpiricalDistribution export BackwardFiniteDifferences export LinearBasisFunctionModel @@ -111,10 +121,14 @@ export BoxBehnken export CentralComposite export CentralFiniteDifferences export CloughPenzien +export DeviationNumber export DoubleLoop export EmpiricalPSD export ExternalModel export SlurmInterface +export ExpectedFeasibility +export ExpectedImprovement +export ExpectedImprovementForGlobalFit export Extractor export FaureSampling export FORM @@ -122,9 +136,11 @@ export ForwardFiniteDifferences export FractionalFactorial export FullFactorial export GaussianMixtureModel +export GaussianProcess export GaussQuadrature export HaltonSampling export HermiteBasis +export IdentityTransformChoice export ImportanceSampling export Interval export IntervalVariable @@ -141,6 +157,9 @@ export LineSampling export SingleComponentMetropolisHastings export MaximumAPosterioriBayesian export MaximumLikelihoodBayesian +export MaximumLikelihoodEstimation +export MaximinDistance +export MaximumVariance export Model export MonomialBasis export MonteCarlo @@ -152,6 +171,7 @@ export PolynomialChaosExpansion export PolyharmonicRadialBasis export PolyharmonicSpline export ProbabilityBox +export ProbabilityOfImprovement export RadialBasedImportanceSampling export GaussianRadialBasis export RandomVariable @@ -161,6 +181,7 @@ export ShinozukaDeodatis export SobolSampling export Solver export SpectralRepresentation +export StandardNormalTransformChoice export StochasticProcessModel export SubSetInfinity export SubSetInfinityAdaptive @@ -170,7 +191,9 @@ export TransportMap export TransportMapFromSamples export TransportMapBayesian export TwoLevelFactorial -export UQTargetDensity +export UnitRangeTransformChoice +export UpperConfidenceBound +export ZScoreTransformChoice # Methods export bayesianupdating @@ -190,6 +213,7 @@ export mapfromdensity export mapfromsamples export mean export multivariate_indices +export optimize_hyperparameters export pdf export periodogram export polynomialchaos @@ -206,6 +230,7 @@ export to_physical_space! export to_standard_normal_space export to_standard_normal_space! export variancediagnostic +export with_gaussian_noise include("util/binning.jl") include("util/fourier-transform.jl") @@ -239,7 +264,13 @@ include("models/model.jl") include("models/imprecise/propagation.jl") include("models/polyharmonicspline.jl") include("models/responsesurface.jl") -include("models//slicingmodel.jl") +include("models/slicingmodel.jl") +include("models/gp/standardization.jl") +include("models/gp/parameterization.jl") +include("models/gp/hyperparametertuning.jl") +include("models/gp/gaussianprocess.jl") +include("models/gp/adaptivegaussianprocess.jl") +include("models/gp/gp_acquisitionfunction.jl") include("models/ipm.jl") include("models/models.jl") diff --git a/src/inputs/inputs.jl b/src/inputs/inputs.jl index 1a1f657b9..753290156 100644 --- a/src/inputs/inputs.jl +++ b/src/inputs/inputs.jl @@ -37,6 +37,8 @@ function names(inputs::Vector{<:UQInput}) return _names end +names(input::UQInput) = only(names([input])) # need this to get the name of a single input in gps + function count_rvs(inputs::Vector{<:UQInput}) random_inputs = filter(i -> isa(i, RandomUQInput) || isa(i, ProbabilityBox), inputs) return mapreduce(dimensions, +, random_inputs) diff --git a/src/models/gp/adaptivegaussianprocess.jl b/src/models/gp/adaptivegaussianprocess.jl new file mode 100644 index 000000000..2f9e643ed --- /dev/null +++ b/src/models/gp/adaptivegaussianprocess.jl @@ -0,0 +1,364 @@ +""" + AdaptiveGaussianProcess( + gp::GP, + input, + model, + output, + acquisition_function, + n_added_points, + n_design_points = 10, + experimental_design = LatinHypercubeSampling(n_design_points); + kwargs... + ) + +Fit a Gaussian-process surrogate from an initial experimental design, then +adaptively enrich its training data with points selected by +`acquisition_function`. + +At each adaptive iteration, candidate points are sampled from `input`, the +acquisition function selects one candidate, `model` is evaluated at that +point, and the Gaussian process is refitted using the enlarged training set. + +# Arguments +- `gp`: Prior Gaussian process specifying the mean function and kernel. +- `input`: Input variable or variables used for the initial design and + candidate sampling. +- `model`: Expensive model evaluated at selected adaptive points. +- `output`: Name of the model output approximated by the surrogate. +- `acquisition_function`: Learning function used to rank candidate points. +- `n_added_points`: Number of adaptively selected training points. +- `n_design_points`: Number of points in the initial experimental design. +- `experimental_design`: Sampling/design method for the initial training data. + +# Keyword Arguments +- `input_transform`: Transformation applied to GP input features. +- `output_transform`: Transformation applied to the GP response. +- `σ²`: Observation-noise variance. +- `learn_noise`: Whether to infer the observation-noise variance. +- `learn_hyperparameters`: Whether to optimize GP hyperparameters on each fit. +- `candidate_sampling`: Monte Carlo sampling method used for acquisition candidates. +- `optimizer`: Hyperparameter-optimization strategy. + +# Examples +```jldoctest +julia> x = RandomVariable(Uniform(-2, 2), :x); + +julia> model = Model(df -> sin.(df.x), :y); + +julia> prior = GP(ZeroMean(), SqExponentialKernel()); + +julia> surrogate = AdaptiveGaussianProcess( + prior, + x, + model, + :y, + MaximumVariance(), + 5, + ); +``` +""" +function AdaptiveGaussianProcess( + gp::GP, + input::Union{UQInput, Vector{<:UQInput}}, + model::Union{UQModel, Vector{<:UQModel}}, + output::Symbol, + acquisition_function::AbstractGaussianProcessAcquisitionFunction, + n_added_points::Int, + n_design_points::Int = 10, + experimental_design::Union{AbstractMonteCarlo, AbstractDesignOfExperiments} = LatinHypercubeSampling( + n_design_points + ); + input_transform::AbstractTransformChoice = IdentityTransformChoice(), + output_transform::AbstractTransformChoice = IdentityTransformChoice(), + σ²::Float64 = 1.0e-10, + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + candidate_sampling::AbstractMonteCarlo = MonteCarlo(100_000), + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation( + Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false) + ), + ) + gp_model = GaussianProcess( + gp, + input, + model, + output; + experimental_design = experimental_design, + input_transform = input_transform, + output_transform = output_transform, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters, + optimizer = optimizer, + ) + + return AdaptiveGaussianProcess( + gp_model, + input, + model, + acquisition_function, + n_added_points, + candidate_sampling = candidate_sampling, + optimizer = optimizer, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters + ) +end + +""" + AdaptiveGaussianProcess( + input, + model, + output, + acquisition_function, + n_added_points, + n_design_points = 10, + experimental_design = LatinHypercubeSampling(n_design_points); + mean_fct = ZeroMean(), + kernel = SqExponentialKernel(), + kwargs... + ) + +Construct and adapt a Gaussian-process surrogate using a zero-mean, +squared-exponential prior by default. Use `mean_fct` and `kernel` to choose a +different prior. +""" +function AdaptiveGaussianProcess( + input::Union{UQInput, Vector{<:UQInput}}, + model::Union{UQModel, Vector{<:UQModel}}, + output::Symbol, + acquisition_function::AbstractGaussianProcessAcquisitionFunction, + n_added_points::Int, + n_design_points::Int = 10, + experimental_design::Union{AbstractMonteCarlo, AbstractDesignOfExperiments} = LatinHypercubeSampling( + n_design_points + ); + mean_fct::AbstractGPs.MeanFunction = ZeroMean(), + kernel::Kernel = SqExponentialKernel(), + input_transform::AbstractTransformChoice = IdentityTransformChoice(), + output_transform::AbstractTransformChoice = IdentityTransformChoice(), + σ²::Float64 = 1.0e-10, + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + candidate_sampling::AbstractMonteCarlo = MonteCarlo(100_000), + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation( + Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false) + ), + ) + + return AdaptiveGaussianProcess( + GP(mean_fct, kernel), + input, + model, + output, + acquisition_function, + n_added_points, + n_design_points, + experimental_design; + input_transform = input_transform, + output_transform = output_transform, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters, + candidate_sampling = candidate_sampling, + optimizer = optimizer, + ) +end + +""" + AdaptiveGaussianProcess( + gp_model, + input, + model, + acquisition_function, + n_added_points; + kwargs... + ) + +Adapt an already-fitted `GaussianProcess` by evaluating `model` at +`n_added_points` selected from candidates sampled from `input`. +""" +function AdaptiveGaussianProcess( + gp_model::GaussianProcess, + input::Union{UQInput, Vector{<:UQInput}}, + model::Union{UQModel, Vector{<:UQModel}}, + acquisition_function::AbstractGaussianProcessAcquisitionFunction, + n_added_points::Int; + candidate_sampling::AbstractMonteCarlo = MonteCarlo(100_000), + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation( + Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false) + ), + σ²::Float64 = gp_model.σ², + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + ) + for i in 1:n_added_points + candidates = sample(input, candidate_sampling) + next_point = _find_next_point(gp_model, candidates, acquisition_function) + + evaluate!(model, next_point) + gp_model = _refit_gp( + gp_model, next_point, optimizer, σ², learn_noise, learn_hyperparameters + ) + + @debug "added point" iteration = i point = NamedTuple(next_point[1, :]) + end + + return gp_model +end + +""" + AdaptiveGaussianProcess( + gp, + data, + input, + model, + output, + acquisition_function, + n_added_points; + kwargs... + ) + +Fit the initial surrogate from `data` with the specified prior `gp`, then +adaptively add points selected from candidates sampled from `input`. The +`model` supplies the expensive response at each selected point. +""" +function AdaptiveGaussianProcess( + gp::GP, + data::DataFrame, + input::Union{UQInput, Vector{<:UQInput}}, + model::Union{UQModel, Vector{<:UQModel}}, + output::Symbol, + acquisition_function::AbstractGaussianProcessAcquisitionFunction, + n_added_points::Int; + input_transform::AbstractTransformChoice = IdentityTransformChoice(), + output_transform::AbstractTransformChoice = IdentityTransformChoice(), + σ²::Float64 = 1.0e-10, + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + candidate_sampling::AbstractMonteCarlo = MonteCarlo(100_000), + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation( + Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false) + ), + ) + gp_model = GaussianProcess( + gp, + data, + output; + input_transform = input_transform, + output_transform = output_transform, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters, + optimizer = optimizer, + ) + + return AdaptiveGaussianProcess( + gp_model, + input, + model, + acquisition_function, + n_added_points; + candidate_sampling = candidate_sampling, + optimizer = optimizer, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters, + ) +end + +""" + AdaptiveGaussianProcess( + data, + input, + model, + output, + acquisition_function, + n_added_points; + mean_fct = ZeroMean(), + kernel = SqExponentialKernel(), + kwargs... + ) + +Fit the initial surrogate from `data`, then adaptively add points selected +from candidates sampled from `input`. Use `mean_fct` and `kernel` to specify +the initial GP prior; `model` is required to evaluate newly selected points. +""" +function AdaptiveGaussianProcess( + data::DataFrame, + input::Union{UQInput, Vector{<:UQInput}}, + model::Union{UQModel, Vector{<:UQModel}}, + output::Symbol, + acquisition_function::AbstractGaussianProcessAcquisitionFunction, + n_added_points::Int; + mean_fct::AbstractGPs.MeanFunction = ZeroMean(), + kernel::Kernel = SqExponentialKernel(), + input_transform::AbstractTransformChoice = IdentityTransformChoice(), + output_transform::AbstractTransformChoice = IdentityTransformChoice(), + σ²::Float64 = 1.0e-10, + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + candidate_sampling::AbstractMonteCarlo = MonteCarlo(100_000), + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation( + Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false) + ), + ) + + return AdaptiveGaussianProcess( + GP(mean_fct, kernel), + data, + input, + model, + output, + acquisition_function, + n_added_points; + input_transform = input_transform, + output_transform = output_transform, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters, + candidate_sampling = candidate_sampling, + optimizer = optimizer, + ) +end + +function _refit_gp( + gp::GaussianProcess, + new_data::DataFrame, + optimizer::AbstractHyperparameterOptimization, + σ²::Float64, + learn_noise::Bool, + learn_hyperparameters::Bool, + ) + # fit Gaussian process to new data + σ² = check_gp_input(σ², learn_noise) + + # Only add point if it is not in df + unique_new_data = antijoin(new_data, gp.training_data; on = names(new_data)) + append!(gp.training_data, unique_new_data) + + # transform data + x = transform(gp.training_data, gp.input_transformer) + y = transform(gp.training_data, gp.output_transformer) + + if learn_hyperparameters + _gp = optimize_hyperparameters( + PriorGP(gp.posterior.prior, σ², learn_noise), x, y, optimizer + ) + σ² = _gp.σ² + + posterior_gp = posterior(_gp(x), y) + else + posterior_gp = posterior(gp.posterior.prior(x, σ²), y) + end + + return GaussianProcess( + posterior_gp, + gp.output, + σ², + gp.input_transformer, + gp.output_transformer, + gp.training_data, + ) +end diff --git a/src/models/gp/gaussianprocess.jl b/src/models/gp/gaussianprocess.jl new file mode 100644 index 000000000..fe90f420e --- /dev/null +++ b/src/models/gp/gaussianprocess.jl @@ -0,0 +1,424 @@ +struct GaussianProcess <: UQModel + posterior::AbstractGPs.PosteriorGP + output::Symbol + σ²::Float64 + input_transformer::GaussianProcessInputTransformer + output_transformer::GaussianProcessOutputTransformer + training_data::DataFrame +end + +function Base.show(io::IO, gp::GaussianProcess) + print(io, "GaussianProcess(") + print(io, "mean=$(gp.posterior.prior.mean), ") + print(io, "kernel=$(gp.posterior.prior.kernel), ") + print(io, "input=$(gp.input_transformer.input), ") + print(io, "output=$(gp.output), ") + print(io, "n_datapoints=$(size(gp.training_data, 1))") + print(io, ")") + return nothing +end + +function Base.show(io::IO, ::MIME"text/plain", gp::GaussianProcess) + println(io, "GaussianProcess") + println(io, " mean: $(gp.posterior.prior.mean)") + println(io, " kernel: $(gp.posterior.prior.kernel)") + println(io, " input: $(gp.input_transformer.input)") + println(io, " output: $(gp.output)") + print(io, " n_datapoints: $(size(gp.training_data, 1))") + return nothing +end + +# function to check the inputs to a GaussianProcess constructor +function check_gp_input(σ²::Float64, learn_noise::Bool) + # check if σ² is ≥0, not using @assert because apparently it can be turned off and shouldn't be used for function input checking (https://discourse.julialang.org/t/efficient-use-of-test-or-assert/75895/4) + if σ² < 0.0 + throw(DomainError(σ², "σ² < 0")) + end + + # σ² should be >0, otherwise the parameterization throws an error + if learn_noise && σ² < eps() + σ² = 1.0e-5 + @warn "learn_noise was set but σ² is too small, setting σ² = $(σ²)" + end + + if !learn_noise && σ² < eps() + @warn "using small σ² < eps() might lead to numerical instabilities" + end + return σ² +end + +""" + GaussianProcess(data::DataFrame, output::Symbol; kwargs...) + +Constructs a `GaussianProcess` model with the specified data and output variable. + +# Arguments +- `data`: A `DataFrame` containing the input and output data. +- `output`: The output variable for the Gaussian process. + +# Keyword Arguments +- `mean_fct`: The mean function for the Gaussian process. Defaults to `ZeroMean`. +- `kernel`: The kernel for the Gaussian process. Defaults to `SqExponentialKernel`. +- `input_transform`: The transformation to apply to the input variables. Defaults to `IdentityTransformChoice`. +- `output_transform`: The transformation to apply to the output variables. Defaults to `IdentityTransformChoice`. +- `σ²`: The noise variance. Defaults to 1.0e-10. +- `learn_noise`: Whether to learn the noise variance. Defaults to `false`. +- `learn_hyperparameters`: Whether to learn the hyperparameters. Defaults to `true`. +- `optimizer`: The optimization algorithm used to learn the hyperparameters. Defaults to `MaximumLikelihoodEstimation(Optim.LBFGS(), Optim.Options(; iterations=100, show_trace=false))`. + +# Examples +```jldoctest +julia> mean_fct = ConstMean(0.0); + +julia> kernel = SqExponentialKernel(); + +julia> data = DataFrame(x = 1:10, y = [1, 4, 10, 15, 24, 37, 50, 62, 80, 101]); + +julia> gp_model = GaussianProcess(data, :y; mean_fct = mean_fct, kernel = kernel, σ² = 1.0e-3); +``` +""" +function GaussianProcess( + data::DataFrame, + output::Symbol; + mean_fct::AbstractGPs.MeanFunction = ZeroMean(), + kernel::Kernel = SqExponentialKernel(), + input_transform::AbstractTransformChoice = IdentityTransformChoice(), + output_transform::AbstractTransformChoice = IdentityTransformChoice(), + σ²::Float64 = 1.0e-10, + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation(Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false)) + ) + + gp = GP(mean_fct, kernel) + + return GaussianProcess( + gp, data, output; + input_transform = input_transform, + output_transform = output_transform, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters, + optimizer = optimizer + ) +end + +""" + GaussianProcess( + gp::GP, + data::DataFrame, + output::Symbol; + kwargs... + ) + +Constructs a Gaussian process model for the given data and output variable using a pre-defined Gaussian process. + +# Arguments +- `gp`: A Gaussian process object, typically from `AbstractGPs`, defining the kernel and mean. +- `data`: A `DataFrame` containing the input and output data. +- `output`: The name of the output (as a `Symbol`) to be modeled as the response variable. + +# Keyword Arguments +- `input_transform`: Choice of transformation that is applied to input features before fitting. + Defaults to [`IdentityTransformChoice()`](@ref). +- `output_transform`: Choice of transformation that is applied to output data before fitting. + Defaults to [`IdentityTransformChoice()`](@ref). +- `σ²`: The noise variance. Defaults to 0.0. +- `learn_noise`: Whether to learn the noise variance. Defaults to false. +- `learn_hyperparameters`: Whether to learn the hyperparameters. Defaults to true. +- `optimizer`: The optimizer for hyperparameter optimization. Defaults to `MaximumLikelihoodEstimation` with `Optim.LBFGS()` and `Optim.Options(; iterations=100, show_trace=false)`. + +# Examples +```jldoctest +julia> gp = GP(0.0, SqExponentialKernel()); + +julia> data = DataFrame(x = 1:10, y = [1, 4, 10, 15, 24, 37, 50, 62, 80, 101]); + +julia> gp_model = GaussianProcess(gp, data, :y); +``` +""" +function GaussianProcess( + gp::GP, + data::DataFrame, + output::Symbol; + input_transform::AbstractTransformChoice = IdentityTransformChoice(), + output_transform::AbstractTransformChoice = IdentityTransformChoice(), + σ²::Float64 = 1.0e-10, + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation(Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false)) + ) + # force learn_noise to false if learn_hyperparameters is false + if !learn_hyperparameters && learn_noise + @warn "learn_hyperparameters is false, setting learn_noise to false" + learn_noise = false + end + σ² = check_gp_input(σ², learn_noise) + + input = propertynames(data[:, Not(output)]) # Is this always the case? + + # build in- and output transforms + input_transformer = fit_input_transform(data, input, input_transform) + output_transformer = fit_output_transform(data, output, output_transform) + + # transform data + x = transform(data, input_transformer) + y = transform(data, output_transformer) + + posterior_gp = nothing + + # optimize hyperparameters + if learn_hyperparameters + _gp = optimize_hyperparameters(PriorGP(gp, σ², learn_noise), x, y, optimizer) + σ² = _gp.σ² + # _gp is a PriorGP object, calling it directly involves the noise, so no need to add σ² again + posterior_gp = posterior(_gp(x), y) + else + # gp has to be called with noise since it is an AbstractGPs.GP object, not a PriorGP object + posterior_gp = posterior(gp(x, σ²), y) + end + + return GaussianProcess( + posterior_gp, + output, + σ², + input_transformer, + output_transformer, + data + ) +end + +""" + GaussianProcess(input::Union{UQInput, Vector{<:UQInput}}, model::Union{UQModel, Vector{<:UQModel}}, output::Symbol; kwargs...) + +Constructs a `GaussianProcess` model with the specified input, model, and output. + +# Arguments +- `input`: The input variable(s) for the Gaussian process. +- `model`: The model(s) to be used for the Gaussian process. +- `output`: The output variable for the Gaussian process. + +# Keyword Arguments +- `n_design_points`: Number of design points to sample from the input space. Defaults to 10. +- `experimental_design`: The strategy utilized for sampling the input variables. Defaults to `LatinHypercubeSampling`. +- `mean_fct`: The mean function for the Gaussian process. Defaults to `ZeroMean`. +- `kernel`: The kernel for the Gaussian process. Defaults to `SqExponentialKernel`. +- `input_transform`: The transformation to apply to the input variables. Defaults to `IdentityTransformChoice`. +- `output_transform`: The transformation to apply to the output variables. Defaults to `IdentityTransformChoice`. +- `σ²`: The noise variance. Defaults to 0.0. +- `learn_noise`: Whether to learn the noise variance. Defaults to `false`. +- `learn_hyperparameters`: Whether to learn the hyperparameters. Defaults to `true`. +- `optimizer`: The optimization algorithm used to learn the hyperparameters. Defaults to `MaximumLikelihoodEstimation(Optim.LBFGS(), Optim.Options(; iterations=100, show_trace=false))`. + +# Examples +```jldoctest +julia> begin # hide + mean_fct = ConstMean(0.0) + kernel = SqExponentialKernel() + x = RandomVariable(Uniform(0, 5), :x) + model = Model(df -> sin.(df.x), :y) + design = LatinHypercubeSampling(10) + gp_model = GaussianProcess(x, model, :y; experimental_design = design, mean_fct = mean_fct, kernel = kernel) + nothing # hide + end # hide +``` +""" +function GaussianProcess( + input::Union{UQInput, Vector{<:UQInput}}, + model::Union{UQModel, Vector{<:UQModel}}, + output::Symbol; + n_design_points::Int = 10, + experimental_design::Union{AbstractMonteCarlo, AbstractDesignOfExperiments} = LatinHypercubeSampling(n_design_points), + mean_fct::AbstractGPs.MeanFunction = ZeroMean(), + kernel::Kernel = SqExponentialKernel(), + input_transform::AbstractTransformChoice = IdentityTransformChoice(), + output_transform::AbstractTransformChoice = IdentityTransformChoice(), + σ²::Float64 = 1.0e-10, + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation(Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false)) + ) + + gp = GP(mean_fct, kernel) + return GaussianProcess( + gp, input, model, output; + experimental_design = experimental_design, + input_transform = input_transform, + output_transform = output_transform, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters, + optimizer = optimizer + ) + +end + +""" + GaussianProcess( + gp::GP, + input::Vector{<:UQInput}, + model::Union{UQModel, Vector{<:UQModel}}, + output::Symbol; + kwargs... + ) + +Constructs a Gaussian process model for the given input and model. Evaluates the model using specified experimental design. + +# Arguments +- `gp`: A Gaussian process object, typically from `AbstractGPs`, defining the kernel and mean. +- `input`: Single input or vector of inputs. The Gaussian process will only consider inputs of type [`RandomVariable`](@ref) as input features. +- `model`: Single model or vector of models of supertype [`UQModel`](@ref) that the Gaussian process is supposed to model. +- `output`: The name of the output (as a `Symbol`) to be modeled as the response variable. + +# Keyword Arguments +- `n_design_points`: Number of design points to sample from the input space. Defaults to 10. +- `experimental_design`: The strategy utilized for sampling the input variables. +- `input_transform`: Choice of transformation that is applied to input features before fitting. + Defaults to [`IdentityTransformChoice()`](@ref). +- `output_transform`: Choice of transformation that is applied to output data before fitting. + Defaults to [`IdentityTransformChoice()`](@ref). +- `σ²`: The noise variance. Defaults to 0.0. +- `learn_noise`: Whether to learn the noise variance. Defaults to `false`. +- `learn_hyperparameters`: Whether to learn the hyperparameters. Defaults to `true`. +- `optimizer`: The optimization algorithm used to learn the hyperparameters. Defaults to `MaximumLikelihoodEstimation(Optim.LBFGS(), Optim.Options(; iterations=100, show_trace=false))`. + +# Examples +```jldoctest +julia> begin # hide + gp = GP(0.0, SqExponentialKernel()) + x = RandomVariable(Uniform(0, 5), :x) + model = Model(df -> sin.(df.x), :y) + design = LatinHypercubeSampling(10) + gp_model = GaussianProcess(gp, x, model, :y; experimental_design = design) + nothing # hide + end # hide +``` +""" +function GaussianProcess( + gp::GP, + input::Vector{<:UQInput}, + model::Union{UQModel, Vector{<:UQModel}}, + output::Symbol; + n_design_points::Int = 10, + experimental_design::Union{AbstractMonteCarlo, AbstractDesignOfExperiments} = LatinHypercubeSampling(n_design_points), + input_transform::AbstractTransformChoice = IdentityTransformChoice(), + output_transform::AbstractTransformChoice = IdentityTransformChoice(), + σ²::Float64 = 1.0e-10, + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation(Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false)) + ) + # build DataFrame + data = sample(input, experimental_design) + evaluate!(model, data) + + # Repeated deterministic input will break the GP kernel + random_input = names(filter(i -> isa(i, RandomVariable), input)) + + return GaussianProcess( + gp, data[!, [random_input..., output]], output; + input_transform = input_transform, + output_transform = output_transform, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters, + optimizer = optimizer + ) +end + +# Helper constructor to wrap `input` into a Vector +function GaussianProcess( + gp::GP, + input::UQInput, + model::Union{UQModel, Vector{<:UQModel}}, + output::Symbol; + n_design_points::Int = 10, + experimental_design::Union{AbstractMonteCarlo, AbstractDesignOfExperiments} = LatinHypercubeSampling(n_design_points), + input_transform::AbstractTransformChoice = IdentityTransformChoice(), + output_transform::AbstractTransformChoice = IdentityTransformChoice(), + σ²::Float64 = 1.0e-10, + learn_noise::Bool = false, + learn_hyperparameters::Bool = true, + optimizer::AbstractHyperparameterOptimization = MaximumLikelihoodEstimation(Optim.LBFGS(), Optim.Options(; iterations = 100, show_trace = false)) + ) + return GaussianProcess( + gp, [input], model, output; + experimental_design = experimental_design, + input_transform = input_transform, + output_transform = output_transform, + σ² = σ², + learn_noise = learn_noise, + learn_hyperparameters = learn_hyperparameters, + optimizer = optimizer + ) +end + +""" + evaluate!(gp::GaussianProcess, data::DataFrame; mode::Symbol = :mean, n_samples::Int = 1) + +Evaluates a fitted [`GaussianProcess`](@ref) model at the specified input locations. + +# Arguments +- `gp`: Trained Gaussian process model to be evaluated. +- `data`: A `DataFrame` containing the input locations at which predictions are computed. + +# Keyword Arguments +- `mode`: A `Symbol` specifying the type of output to return. + Supported options are: + - `:mean` - predictive mean (default) + - `:var` - predictive variance + - `:mean_and_var` - both mean and variance + - `:sample` - random samples from the predictive distribution +- `n_samples`: Number of samples to draw when `mode = :sample`. Ignored otherwise. + (Note: Sampling can be unstable when input locations are very close together, leading to numerical issues in the covariance matrix.) + +# Examples +```jldoctest +julia> gp = GP(0.0, SqExponentialKernel()); + +julia> data = DataFrame(x = 1:10, y = [1, 4, 10, 15, 24, 37, 50, 62, 80, 101]); + +julia> gp_model = GaussianProcess(gp, data, :y; σ² = 1.0e-3); + +julia> df = DataFrame(x = [0.5, 1.5, 2.5, 5.5, 8.5]); + +julia> evaluate!(gp_model, df; mode = :mean_and_var); +``` +""" +function evaluate!( + gp::GaussianProcess, + data::DataFrame; + mode::Symbol = :mean, + n_samples::Int = 1 + ) + x = transform(data, gp.input_transformer) + finite_projection = gp.posterior(x, gp.σ²) + + if mode === :mean + μ = mean(finite_projection) + col = Symbol(string(gp.output, "_mean")) + data[!, col] = inverse_transform(μ, gp.output_transformer) + elseif mode === :var + σ² = var(finite_projection) + col = Symbol(string(gp.output, "_var")) + data[!, col] = variance_inverse_transform(σ², gp.output_transformer) + elseif mode === :mean_and_var + μ = mean(finite_projection) + σ² = var(finite_projection) + col_mean = Symbol(string(gp.output, "_mean")) + col_var = Symbol(string(gp.output, "_var")) + data[!, col_mean] = inverse_transform(μ, gp.output_transformer) + data[!, col_var] = variance_inverse_transform(σ², gp.output_transformer) + elseif mode === :sample + samples = rand(finite_projection, n_samples) + cols = [Symbol(string(gp.output, "_sample_", i)) for i in 1:n_samples] + foreach( + (col, sample) -> data[!, col] = inverse_transform(sample, gp.output_transformer), + cols, eachcol(samples) + ) + else + throw(ArgumentError("Unknown `GaussianProcess` evaluation mode: $mode")) + end + + return nothing +end diff --git a/src/models/gp/gp_acquisitionfunction.jl b/src/models/gp/gp_acquisitionfunction.jl new file mode 100644 index 000000000..569f59864 --- /dev/null +++ b/src/models/gp/gp_acquisitionfunction.jl @@ -0,0 +1,319 @@ +#=== +The type `AbstractGaussianProcessAcquisitionFunction` is used to later be able to implement +acquisition functions for other surrogates, e.g., PCE. + +Some common acquisition (learning) functions for adaptive GPs for different kind of applications + +### General Active Learning +Goal: Get good fit of the GP (globally) +- `MaximumVariance()` simply adds a new point based on the maximum variance +Some other learning functions are implemented based on Fuhg et al. (2021) and they provide a +tradeoff between exploitation (refine local extrema) and exploration (reduce global variance). +Specifically, there are the more exploration-based functions: +- `MaximinDistance()` +- `ExpectedImprovementForGlobalFit` +There are more learning functions that can be implemented (https://github.com/FuhgJan/StateOfTheArtAdaptiveSampling). + +### Bayesian Optimization +Goal: Refine the global minimum of the function +- `ExpectedImprovement()` +- `ProbabilityOfImprovement()` +- `UpperConfidenceBound()` + +### Reliability Analysis +Goal: Refine limit state surface, `g(x) = 0`, +- `DeviationNumber()` (`U`-function from AK-MCS) +- `ExpectedFeasibility()` (`EFF` from EGRA) + +### References +Fuhg, J. N., Fau, A., & Nackenhorst, U. (2021). State-of-the-Art and Comparative Review +of Adaptive Sampling Methods for Kriging. Archives of Computational Methods in Engineering, +28(4), 2689–2747. https://doi.org/10.1007/s11831-020-09474-6 +===# + +""" + MaximumVariance() + +Selects the candidate with the largest posterior variance. + +### Reference +Sacks, J., Welch, W. J., Mitchell, T. J., & Wynn, H. P. (1989). Design and +Analysis of Computer Experiments. Statistical Science, 4(4), 409–423. +https://doi.org/10.1214/ss/1177012413 +""" +struct MaximumVariance <: AbstractGaussianProcessAcquisitionFunction end + +function _find_next_point(gp::GaussianProcess, candidates::DataFrame, ::MaximumVariance) + input = propertynames(candidates) + candidates = copy(candidates) + _, var_col = _mean_var_cols(gp) + + evaluate!(gp, candidates; mode = :var) + ind = argmax(candidates[:, var_col]) + return candidates[[ind], Cols(input)] +end + +""" + ExpectedImprovement(ξ = 0.0) + +Selects the candidate maximizing expected improvement over the current best +(minimum) observed training output. `ξ` controls the exploration/exploitation +trade-off (higher `ξ` favors exploration). + +### Reference +Jones, D. R., Schonlau, M., & Welch, W. J. (1998). Efficient Global +Optimization of Expensive Black-Box Functions. Journal of Global +Optimization, 13(4), 455–492. https://doi.org/10.1023/A:1008306431147 +""" +Base.@kwdef struct ExpectedImprovement <: AbstractGaussianProcessAcquisitionFunction + ξ::Float64 = 0.0 +end + +function _find_next_point(gp::GaussianProcess, candidates::DataFrame, ei::ExpectedImprovement) + input = propertynames(candidates) + candidates = copy(candidates) + mean_col, var_col = _mean_var_cols(gp) + + evaluate!(gp, candidates; mode = :mean_and_var) + μ = candidates[:, mean_col] + σ = sqrt.(candidates[:, var_col]) + f_best = minimum(gp.training_data[:, gp.output]) + + improvement = f_best .- μ .- ei.ξ + z = improvement ./ σ + ei_values = improvement .* cdf.(Normal(), z) .+ σ .* pdf.(Normal(), z) + ei_values[σ .<= 0] .= 0.0 + + ind = argmax(ei_values) + return candidates[[ind], Cols(input)] +end + +""" + ProbabilityOfImprovement(ξ = 0.0) + +Selects the candidate maximizing the probability of improving over the current +best (minimum) observed training output. + +### Reference +Kushner, H. J. (1964). A New Method of Locating the Maximum Point of an +Arbitrary Multipeak Curve in the Presence of Noise. Journal of Basic +Engineering, 86(1), 97–106. https://doi.org/10.1115/1.3653121 +""" +Base.@kwdef struct ProbabilityOfImprovement <: AbstractGaussianProcessAcquisitionFunction + ξ::Float64 = 0.0 +end + +function _find_next_point(gp::GaussianProcess, candidates::DataFrame, poi::ProbabilityOfImprovement) + input = propertynames(candidates) + candidates = copy(candidates) + mean_col, var_col = _mean_var_cols(gp) + + evaluate!(gp, candidates; mode = :mean_and_var) + μ = candidates[:, mean_col] + σ = sqrt.(candidates[:, var_col]) + f_best = minimum(gp.training_data[:, gp.output]) + + z = (f_best .- μ .- poi.ξ) ./ σ + pi_values = cdf.(Normal(), z) + pi_values[σ .<= 0] .= 0.0 + + ind = argmax(pi_values) + return candidates[[ind], Cols(input)] +end + + +""" + UpperConfidenceBound(κ = 2.0) + +Selects the candidate minimizing `μ(x) - κ·σ(x)` (lower confidence bound for a +minimization objective). `κ` controls the exploration weight. + +### Reference +Cox, D. D., & John, S. (1992). A Statistical Method for Global Optimization. +Proceedings of the 1992 IEEE International Conference on Systems, Man, and +Cybernetics, 1241–1246. https://doi.org/10.1109/ICSMC.1992.271617 +""" +Base.@kwdef struct UpperConfidenceBound <: AbstractGaussianProcessAcquisitionFunction + κ::Float64 = 2.0 +end + +function _find_next_point(gp::GaussianProcess, candidates::DataFrame, ucb::UpperConfidenceBound) + input = propertynames(candidates) + candidates = copy(candidates) + mean_col, var_col = _mean_var_cols(gp) + + evaluate!(gp, candidates; mode = :mean_and_var) + μ = candidates[:, mean_col] + σ = sqrt.(candidates[:, var_col]) + + lcb = μ .- ucb.κ .* σ + ind = argmin(lcb) + return candidates[[ind], Cols(input)] +end + +""" + DeviationNumber(threshold = 0.0) + +U-function: selects the candidate minimizing `|μ(x) - threshold| / σ(x)`, i.e. the +point closest to the limit state and with the highest uncertainty. + +### Reference +Echard, B., Gayton, N., & Lemaire, M. (2011). AK-MCS: An active learning +reliability method combining Kriging and Monte Carlo Simulation. Structural Safety, 33(2), +145–154. https://doi.org/10.1016/j.strusafe.2011.01.002 +""" +Base.@kwdef struct DeviationNumber <: AbstractGaussianProcessAcquisitionFunction + threshold::Float64 = 0.0 + stopping::Float64 = 2.0 +end + +function _find_next_point(gp::GaussianProcess, candidates::DataFrame, dn::DeviationNumber) + next_point, _ = _find_next_point_stopping(gp, candidates, dn) + return next_point +end + +function _find_next_point_stopping(gp::GaussianProcess, candidates::DataFrame, dn::DeviationNumber) + input = propertynames(candidates) + candidates = copy(candidates) + mean_col, var_col = _mean_var_cols(gp) + + evaluate!(gp, candidates; mode = :mean_and_var) + μ = candidates[:, mean_col] + σ = sqrt.(candidates[:, var_col]) + + u = abs.(μ .- dn.threshold) ./ σ + val, ind = findmin(u) + + return candidates[[ind], Cols(input)], (val >= dn.stopping) +end + +function _mean_var_cols(gp::GaussianProcess) + return Symbol(gp.output, "_mean"), Symbol(gp.output, "_var") +end + +""" + ExpectedFeasibility(threshold = 0.0, epsilon_factor = 2.0) + +Selects the candidate with the largest expected feasibility function (EFF) +near the limit state `G(x) = threshold`. The feasibility half-width is set +per candidate as `epsilon_factor * σ(x)`. + +### Reference +Bichon, B. J., Eldred, M. S., Swiler, L. P., Mahadevan, S., & McFarland, +J. M. (2008). Efficient Global Reliability Analysis for Nonlinear Implicit +Performance Functions. AIAA Journal, 46(10), 2459-2468. +https://doi.org/10.2514/1.34321 +""" +Base.@kwdef struct ExpectedFeasibility <: AbstractGaussianProcessAcquisitionFunction + threshold::Float64 = 0.0 + epsilon_factor::Float64 = 2.0 + stopping::Float64 = 0.001 +end + +function _find_next_point( + gp::GaussianProcess, + candidates::DataFrame, + eff::ExpectedFeasibility, + ) + next_point, _ = _find_next_point_stopping(gp, candidates, eff) + return next_point +end + +function _find_next_point_stopping( + gp::GaussianProcess, + candidates::DataFrame, + eff::ExpectedFeasibility, + ) + input = propertynames(candidates) + candidates = copy(candidates) + mean_col, var_col = _mean_var_cols(gp) + + evaluate!(gp, candidates; mode = :mean_and_var) + μ = candidates[:, mean_col] + σ = sqrt.(candidates[:, var_col]) + + eff_values = zeros(length(μ)) + positive_σ = σ .> 0 + + μ_active = μ[positive_σ] + σ_active = σ[positive_σ] + ε = eff.epsilon_factor .* σ_active + + z = (eff.threshold .- μ_active) ./ σ_active + z_lower = (eff.threshold .- ε .- μ_active) ./ σ_active + z_upper = (eff.threshold .+ ε .- μ_active) ./ σ_active + + eff_values[positive_σ] = + (μ_active .- eff.threshold) .* ( + 2 .* cdf.(Normal(), z) .- + cdf.(Normal(), z_lower) .- + cdf.(Normal(), z_upper) + ) .- + σ_active .* ( + 2 .* pdf.(Normal(), z) .- + pdf.(Normal(), z_lower) .- + pdf.(Normal(), z_upper) + ) .+ + ε .* ( + cdf.(Normal(), z_upper) .- + cdf.(Normal(), z_lower) + ) + + val, ind = findmax(eff_values) + + return candidates[[ind], Cols(input)], (val <= eff.stopping) +end + +""" + MaximinDistance() + +Space-filling score that selects the candidate farthest from every +existing training point (maximizes the nearest-neighbor distance). + +### Reference +Johnson, M. E., Moore, L. M., & Ylvisaker, D. (1990). Minimax and Maximin +Distance Designs. Journal of Statistical Planning and Inference, 26(2), +131–148. https://doi.org/10.1016/0378-3758(90)90122-B +""" +struct MaximinDistance <: AbstractGaussianProcessAcquisitionFunction end + +function _find_next_point(gp::GaussianProcess, candidates::DataFrame, ::MaximinDistance) + input = propertynames(candidates) + X = Matrix(gp.training_data[:, input]) + Xc = Matrix(candidates[:, input]) + + distances = [minimum(norm(Xc[i, :] - X[j, :]) for j in axes(X, 1)) for i in axes(Xc, 1)] + ind = argmax(distances) + return candidates[[ind], Cols(input)] +end + +""" + ExpectedImprovementForGlobalFit() + +EIGF score rewards candidates far (in output) from their nearest +training observation, adjusted by local posterior variance. + +### Reference +Lam, C. Q. (2008). Sequential Adaptive Designs in Computer Experiments for +Response Surface Model Fit. PhD dissertation, The Ohio State University. +""" +struct ExpectedImprovementForGlobalFit <: AbstractGaussianProcessAcquisitionFunction end + +function _find_next_point(gp::GaussianProcess, candidates::DataFrame, ::ExpectedImprovementForGlobalFit) + input = propertynames(candidates) + candidates = copy(candidates) + mean_col, var_col = _mean_var_cols(gp) + evaluate!(gp, candidates; mode = :mean_and_var) + + μ = candidates[:, mean_col] + σ² = candidates[:, var_col] + X = Matrix(gp.training_data[:, input]) + Xc = Matrix(candidates[:, input]) + y = gp.training_data[:, gp.output] + + nearest = [argmin([norm(Xc[i, :] - X[j, :]) for j in axes(X, 1)]) for i in axes(Xc, 1)] + eigf = abs2.(μ .- y[nearest]) .+ σ² + + ind = argmax(eigf) + return candidates[[ind], Cols(input)] +end diff --git a/src/models/gp/hyperparametertuning.jl b/src/models/gp/hyperparametertuning.jl new file mode 100644 index 000000000..2a1662a9a --- /dev/null +++ b/src/models/gp/hyperparametertuning.jl @@ -0,0 +1,125 @@ +abstract type AbstractHyperparameterOptimization end + +""" + MaximumLikelihoodEstimation( + optimizer = LBFGS(), + options = Optim.Options(; iterations = 100, show_trace = false), + backend = AutoMooncake(); + restarts = 5, + ) + +Represents a hyperparameter optimization strategy that maximizes the log marginal likelihood +of a Gaussian process model with random restarts of the optimization. + +# Arguments +- `optimizer::Optim.AbstractOptimizer`: chosen optimizer (default: `LBFGS()`) +- `options::Optim.Options`: options for the optimizer (default: `Optim.Options(; iterations = 100, show_trace = false)`) +- `backend::AbstractADType`: automatic differentiation backend (default: `AutoMooncake()`; can be `nothing` when using gradient-free optimization) +# Keyword Arguments +- `restarts`: Number of additional randomized optimization runs. Defaults to `5`. + +# Note +You can choose from any optimizer and set of options provided by [`Optim.jl`](https://julianlsolvers.github.io/Optim.jl/stable/), +such as `LBFGS()`, `Adam()`, or `ConjugateGradient()`. + +# Examples + +```jldoctest +julia> using Optim + +julia> MaximumLikelihoodEstimation(Optim.Adam(alpha = 0.01), Optim.Options(; iterations = 1000, show_trace = false)) +MaximumLikelihoodEstimation(Adam{Float64, Float64, Flat}(0.01, 0.9, 0.999, 1.0e-8, Flat()), Optim.Options(x_abstol = 0.0, x_reltol = 0.0, f_abstol = 0.0, f_reltol = 0.0, g_abstol = 1.0e-8, outer_x_abstol = 0.0, outer_x_reltol = 0.0, outer_f_abstol = 0.0, outer_f_reltol = 0.0, outer_g_abstol = 1.0e-8, f_calls_limit = 0, g_calls_limit = 0, h_calls_limit = 0, allow_f_increases = true, allow_outer_f_increases = true, successive_f_tol = 1, iterations = 1000, outer_iterations = 1000, store_trace = false, trace_simplex = false, show_trace = false, extended_trace = false, show_warnings = true, show_every = 1, time_limit = NaN, ) +, AutoMooncake(), 5) +``` +""" +struct MaximumLikelihoodEstimation <: AbstractHyperparameterOptimization + optimizer::Optim.AbstractOptimizer + options::Optim.Options + backend::Union{AbstractADType, Nothing} + restarts::Int + + function MaximumLikelihoodEstimation( + optimizer::Optim.AbstractOptimizer = LBFGS(), + options::Optim.Options = Optim.Options(; iterations = 100, show_trace = false), + backend::Union{AbstractADType, Nothing} = AutoMooncake(); + restarts::Int = 5, + ) + restarts >= 0 || throw(ArgumentError("restarts must be nonnegative")) + + if isa(optimizer, Optim.ZerothOrderOptimizer) + return new(optimizer, options, nothing, restarts) + end + + return new(optimizer, options, backend, restarts) + end +end + +objective( + f::PriorGP, + x::Union{RowVecs{<:Real}, Vector{<:Real}}, + y::Vector{<:Real}, + ::MaximumLikelihoodEstimation +) = -logpdf(f(x), y) + +_initializer(θ) = θ .+ 0.5 .* randn(length(θ)) + +function optimize_hyperparameters( + gp::PriorGP, + x::Union{RowVecs{<:Real}, Vector{<:Real}}, + y::Vector{<:Real}, + mle::MaximumLikelihoodEstimation + ) + model, θ₀ = parameterize(gp) + θ₀_flat, unflatten = ParameterHandling.flatten(θ₀) + obj = θ -> objective(model(unflatten(θ)), x, y, mle) + + best_gp = nothing + best_objective = Inf + + for restart in 0:mle.restarts + θ_init = restart == 0 ? copy(θ₀_flat) : _initializer(θ₀_flat) + gp_opt = _optimize_hyperparameters(obj, model, unflatten, θ_init, mle) + value = objective(gp_opt, x, y, mle) + + if value < best_objective + best_objective = value + best_gp = gp_opt + end + end + + return best_gp +end + +function _optimize_hyperparameters( + obj::Function, + model, + unflatten::Function, + θ_init::AbstractVector, + mle::MaximumLikelihoodEstimation, + ) + if isa(mle.optimizer, Optim.FirstOrderOptimizer) || + isa(mle.optimizer, Optim.SecondOrderOptimizer) + + prep = DifferentiationInterface.prepare_gradient(obj, mle.backend, θ_init) + + function fg!(F, G, θ) + value, gradient = DifferentiationInterface.value_and_gradient( + obj, prep, mle.backend, θ + ) + G !== nothing && (G .= gradient) + return value + end + + result = optimize(NLSolversBase.only_fg!(fg!), θ_init, mle.optimizer, mle.options) + + return model(unflatten(result.minimizer)) + + elseif isa(mle.optimizer, Optim.ZerothOrderOptimizer) + # Gradient-free optimizer + result = optimize(obj, θ_init, mle.optimizer, mle.options) + + return model(unflatten(result.minimizer)) + else + error("Optimizer of type $(typeof(mle.optimizer)) not supported.") + end +end diff --git a/src/models/gp/parameterization.jl b/src/models/gp/parameterization.jl new file mode 100644 index 000000000..a536bc90c --- /dev/null +++ b/src/models/gp/parameterization.jl @@ -0,0 +1,263 @@ +""" +# Developer Note + +`Parameterized(object)` wraps `object` so it can be called with parameters `θ`. + +`parameterize(object)` returns a parameterized, callable version of the object and its parameters. + ```julia + model, θ = parameterize(obj) + model(θ) # returns a new object with parameters applied + This works for mean functions, kernels, transformations, and Gaussian processes. + +Based on two core functions, this system can extract model parameters for an optimization routine +and apply potentially constrained parameters to the underlying model to compute the optimization objective. +The two core functions are: + + 1. extract_parameters(obj) + + Returns the free parameters of obj wrapped in ParameterHandling containers. + Enforces constraints (e.g., positive or bounded) where applicable. + For composite objects (like a `AbstractGPs.GP`), returns a tuple of componentwise parameter sets. + Returns nothing for objects without trainable parameters. + + 2. apply_parameters(obj, θ) + + Returns a new object of the same type with parameters θ applied. + For hierarchical objects, θ is expected to match the structure returned by extract_parameters. +""" + +struct Parameterized{T} + object::T +end + +function (p::Parameterized)(θ) + return apply_parameters(p.object, ParameterHandling.value(θ)) +end + +parameterize(object) = Parameterized(object), extract_parameters(object) + +extract_parameters(::ZeroMean) = nothing +apply_parameters(m::ZeroMean, θ) = m + +extract_parameters(m::ConstMean) = m.c +apply_parameters(::ConstMean, θ) = ConstMean(θ) + +extract_parameters(::CustomMean) = nothing +apply_parameters(m::CustomMean, θ) = m + +extract_parameters(::ZeroKernel) = nothing +apply_parameters(k::ZeroKernel, _) = k + +extract_parameters(::WhiteKernel) = nothing +apply_parameters(k::WhiteKernel, _) = k + +extract_parameters(::CosineKernel) = nothing +apply_parameters(k::CosineKernel, _) = k + +extract_parameters(::SqExponentialKernel) = nothing +apply_parameters(k::SqExponentialKernel, _) = k + +extract_parameters(::ExponentialKernel) = nothing +apply_parameters(k::ExponentialKernel, _) = k + +extract_parameters(::ExponentiatedKernel) = nothing +apply_parameters(k::ExponentiatedKernel, _) = k + +extract_parameters(::Matern32Kernel) = nothing +apply_parameters(k::Matern32Kernel, _) = k + +extract_parameters(::Matern52Kernel) = nothing +apply_parameters(k::Matern52Kernel, _) = k + +extract_parameters(::Matern72Kernel) = nothing +apply_parameters(k::Matern72Kernel, _) = k + +extract_parameters(::NeuralNetworkKernel) = nothing +apply_parameters(k::NeuralNetworkKernel, _) = k + +extract_parameters(::PiecewisePolynomialKernel) = nothing +apply_parameters(k::PiecewisePolynomialKernel, _) = k + +extract_parameters(::WienerKernel) = nothing +apply_parameters(k::WienerKernel, _) = k + +extract_parameters(::FunctionTransform) = nothing +apply_parameters(t::FunctionTransform, _) = t + +extract_parameters(::SelectTransform) = nothing +apply_parameters(t::SelectTransform, _) = t + +extract_parameters(::IdentityTransform) = nothing +apply_parameters(t::IdentityTransform, _) = t + +function extract_parameters(::IndependentMOKernel) + throw( + ArgumentError("IndependentMOKernel not supported for hyper parameter optimization.") + ) +end +function apply_parameters(::IndependentMOKernel, _) + throw( + ArgumentError("IndependentMOKernel not supported for hyper parameter optimization.") + ) +end +function extract_parameters(::IntrinsicCoregionMOKernel) + throw( + ArgumentError( + "IntrinsicCoregionMOKernel not supported for hyper parameter optimization." + ), + ) +end +function apply_parameters(::IntrinsicCoregionMOKernel, _) + throw( + ArgumentError( + "IntrinsicCoregionMOKernel not supported for hyper parameter optimization." + ), + ) +end +function extract_parameters(::LatentFactorMOKernel) + throw(ArgumentError("LatentFactorMOKernel not supported hyper parameter optimization.")) +end +function apply_parameters(::LatentFactorMOKernel, _) + throw( + ArgumentError( + "LatentFactorMOKernel not supported for hyper parameter optimization." + ), + ) +end +function extract_parameters(::LinearMixingModelKernel) + throw( + ArgumentError("LinearMixingModelKernel not supported hyper parameter optimization.") + ) +end +function apply_parameters(::LinearMixingModelKernel, _) + throw( + ArgumentError( + "LinearMixingModelKernel not supported for hyper parameter optimization." + ), + ) +end +function extract_parameters(::KernelFunctions.NeuralKernelNetwork) + throw(ArgumentError("NeuralKernelNetwork not supported hyper parameter optimization.")) +end +function apply_parameters(::KernelFunctions.NeuralKernelNetwork, _) + throw( + ArgumentError("NeuralKernelNetwork not supported for hyper parameter optimization.") + ) +end +function extract_parameters(::GibbsKernel) + throw(ArgumentError("GibbsKernel not supported hyper parameter optimization.")) +end +function apply_parameters(::GibbsKernel, _) + throw( + ArgumentError("GibbsKernel not supported for hyper parameter optimization.") + ) +end + +extract_parameters(k::ConstantKernel) = ParameterHandling.positive(k.c) +apply_parameters(::ConstantKernel, θ) = ConstantKernel(; c = only(θ)) + +extract_parameters(k::GammaExponentialKernel) = ParameterHandling.bounded(k.γ, 0.0, 2.0) +apply_parameters(::GammaExponentialKernel, θ) = GammaExponentialKernel(; γ = only(θ)) + +extract_parameters(k::FBMKernel) = ParameterHandling.bounded(k.h, 0.0, 1.0) +apply_parameters(::FBMKernel, θ) = FBMKernel(; h = only(θ)) + +extract_parameters(k::MaternKernel) = ParameterHandling.positive(k.ν) +apply_parameters(::MaternKernel, θ) = MaternKernel(; ν = only(θ)) + +extract_parameters(k::PeriodicKernel) = ParameterHandling.positive(k.r) +apply_parameters(::PeriodicKernel, θ) = PeriodicKernel(; r = θ) + +extract_parameters(k::LinearKernel) = ParameterHandling.positive(k.c) +apply_parameters(::LinearKernel, θ) = LinearKernel(; c = only(θ)) + +extract_parameters(k::PolynomialKernel) = ParameterHandling.positive(k.c) +apply_parameters(::PolynomialKernel, θ) = PolynomialKernel(; c = only(θ)) + +extract_parameters(k::RationalKernel) = ParameterHandling.positive(k.α) +apply_parameters(::RationalKernel, θ) = RationalKernel(; α = only(θ)) + +extract_parameters(k::RationalQuadraticKernel) = ParameterHandling.positive(k.α) +apply_parameters(::RationalQuadraticKernel, θ) = RationalQuadraticKernel(; α = only(θ)) + +function extract_parameters(k::GammaRationalKernel) + return (ParameterHandling.positive(k.α), ParameterHandling.bounded(k.γ, 0.0, 2.0)) +end +function apply_parameters(::GammaRationalKernel, θ) + return GammaRationalKernel(; α = only(θ[1]), γ = only(θ[2])) +end + +# kernels (see KernelFunctions.jl src/kernels) +extract_parameters(k::KernelProduct) = map(extract_parameters, k.kernels) +apply_parameters(k::KernelProduct, θ) = KernelProduct(map(apply_parameters, k.kernels, θ)) + +extract_parameters(k::KernelSum) = map(extract_parameters, k.kernels) +apply_parameters(k::KernelSum, θ) = KernelSum(map(apply_parameters, k.kernels, θ)) + +extract_parameters(k::KernelTensorProduct) = map(extract_parameters, k.kernels) +function apply_parameters(k::KernelTensorProduct, θ) + return KernelTensorProduct(map(apply_parameters, k.kernels, θ)) +end + +extract_parameters(k::NormalizedKernel) = extract_parameters(k.kernel) +apply_parameters(k::NormalizedKernel, θ) = NormalizedKernel(apply_parameters(k.kernel, θ)) + +function extract_parameters(k::ScaledKernel) + return (extract_parameters(k.kernel), ParameterHandling.positive(only(k.σ²))) +end +apply_parameters(k::ScaledKernel, θ) = ScaledKernel(apply_parameters(k.kernel, θ[1]), θ[2]) + +function extract_parameters(k::TransformedKernel) + return (extract_parameters(k.kernel), extract_parameters(k.transform)) +end +function apply_parameters(k::TransformedKernel, θ) + return TransformedKernel(apply_parameters(k.kernel, θ[1]), apply_parameters(k.transform, θ[2])) +end + +# transform (see KernelFunctions.jl src/transform) +extract_parameters(t::ARDTransform) = ParameterHandling.positive(t.v) +apply_parameters(::ARDTransform, θ) = ARDTransform(θ) + +extract_parameters(t::ChainTransform) = map(extract_parameters, t.transforms) +function apply_parameters(t::ChainTransform, θ) + return ChainTransform(map(apply_parameters, t.transforms, θ)) +end + +extract_parameters(t::LinearTransform) = t.A +apply_parameters(::LinearTransform, θ) = LinearTransform(θ) + +extract_parameters(t::PeriodicTransform) = ParameterHandling.positive(t.f) +apply_parameters(::PeriodicTransform, θ) = PeriodicTransform(θ) + +extract_parameters(t::ScaleTransform) = ParameterHandling.positive(t.s) +apply_parameters(::ScaleTransform, θ) = ScaleTransform(θ) + +# ---------------- Gaussian Processes ---------------- +extract_parameters(f::GP) = (extract_parameters(f.mean), extract_parameters(f.kernel)) +function apply_parameters(f::GP, θ) + return GP(apply_parameters(f.mean, θ[1]), apply_parameters(f.kernel, θ[2])) +end + + +#Internal struct for hyperparameter optimization. The struct saves the GP that is to be optimized and the noise. +struct PriorGP{T <: GP, Tn <: Real} + gp::T + σ²::Tn + learn_noise::Bool +end + +(gp::PriorGP)(x) = gp.gp(x, gp.σ²) + +extract_parameters(f::PriorGP) = begin + f.learn_noise ? ( + extract_parameters(f.gp), + ParameterHandling.positive(f.σ², exp, 1.0e-6), + ) : + (extract_parameters(f.gp)) +end + +apply_parameters(f::PriorGP, θ) = begin + f.learn_noise ? + PriorGP(apply_parameters(f.gp, θ[1]), θ[2], f.learn_noise) : + PriorGP(apply_parameters(f.gp, θ), f.σ², f.learn_noise) +end diff --git a/src/models/gp/standardization.jl b/src/models/gp/standardization.jl new file mode 100644 index 000000000..3965d10a8 --- /dev/null +++ b/src/models/gp/standardization.jl @@ -0,0 +1,267 @@ +abstract type AbstractTransformChoice end + +""" + IdentityTransformChoice() + +A standardization choice that specifies the application of an identity transform to data. + +Used as an input or output transformation in a [`GaussianProcess`](@ref). +Internally, the `DataStandardizer` constructs the functions required for evaluation. + +# Examples +```jldoctest +julia> id = IdentityTransformChoice() +IdentityTransformChoice() +``` +""" +struct IdentityTransformChoice <: AbstractTransformChoice end + +""" + ZScoreTransformChoice() + +A standardization choice that specifies the application of a Z-score-transformation to data. + +Used as an input or output transformation in a [`GaussianProcess`](@ref). +Internally, the `DataStandardizer` constructs the functions required for evaluation. + +# Examples +```jldoctest +julia> zscore = ZScoreTransformChoice() +ZScoreTransformChoice() +``` +""" +struct ZScoreTransformChoice <: AbstractTransformChoice end + +""" + UnitRangeTransformChoice() + +A standardization choice that specifies the application of a unit range transform to data. + +Used as an input or output transformation in a [`GaussianProcess`](@ref). +Internally, the `DataStandardizer` constructs the functions required for evaluation. + +# Examples +```jldoctest +julia> unitrange = UnitRangeTransformChoice() +UnitRangeTransformChoice() +``` +""" +struct UnitRangeTransformChoice <: AbstractTransformChoice end + +""" + StandardNormalTransformChoice() + +A standardization choice that specifies the application of a standard normal transformation to data. + +Can only be used as an input transformation in a [`GaussianProcess`](@ref) for inputs of type [`RandomVariable`](@ref). +Internally, the `DataStandardizer` constructs the function required for evaluation. + +# Examples +```jldoctest +julia> sns = StandardNormalTransformChoice() +StandardNormalTransformChoice() +``` +""" +struct StandardNormalTransformChoice <: AbstractTransformChoice end + + +# ---------------- Utility ---------------- +to_gp_format(x::Vector) = x +to_gp_format(x::Matrix) = RowVecs(x) +dataframe_to_array(df::DataFrame, name::Symbol) = df[:, name] +dataframe_to_array(df::DataFrame, names::Vector{<:Symbol}) = length(names) == 1 ? x = dataframe_to_array(df, only(names)) : x = Matrix(df[:, names]) + +# Internal transform types for dispatching +struct NoTransform end +struct StandardNormalTransform end + + +# ---------------- Input transformation ---------------- +# # Developer Note +# Gaussian process regression inputs are always transformed to Vector or RowVecs in the multivariate input case +struct GaussianProcessInputTransformer{T} + transform::T + input::Union{Symbol, Vector{<:Symbol}, UQInput, Vector{<:UQInput}} +end + +# Fitting +fit_input_transform( + ::DataFrame, + input::Union{Symbol, Vector{<:Symbol}}, + ::IdentityTransformChoice +) = GaussianProcessInputTransformer(NoTransform(), input) + +function fit_input_transform( + data::DataFrame, + input::Union{Symbol, Vector{<:Symbol}}, + ::ZScoreTransformChoice + ) + transform = fit( + StatsBase.ZScoreTransform, + dataframe_to_array(data, input); + dims = 1 + ) + return GaussianProcessInputTransformer(transform, input) +end + +function fit_input_transform( + data::DataFrame, + input::Union{Symbol, Vector{<:Symbol}}, + ::UnitRangeTransformChoice + ) + transform = fit( + StatsBase.UnitRangeTransform, + dataframe_to_array(data, input); + dims = 1 + ) + return GaussianProcessInputTransformer(transform, input) +end + +fit_input_transform( + ::DataFrame, + input::Union{Symbol, Vector{<:Symbol}}, + ::StandardNormalTransformChoice +) = throw(ArgumentError("Standard normal input transform is only valid for inputs of type UQInput")) + +fit_input_transform( + data::DataFrame, + input::Union{UQInput, Vector{<:UQInput}}, + choice::Union{IdentityTransformChoice, ZScoreTransformChoice, UnitRangeTransformChoice} +) = fit_input_transform(data, names(input), choice) + +fit_input_transform( + ::DataFrame, + input::Union{UQInput, Vector{<:UQInput}}, + ::StandardNormalTransformChoice +) = GaussianProcessInputTransformer(StandardNormalTransform(), input) + +# Transforms +transform( + data::DataFrame, + transformer::GaussianProcessInputTransformer{<:NoTransform} +) = to_gp_format(dataframe_to_array(data, transformer.input)) + +transform( + data::DataFrame, + transformer::Union{GaussianProcessInputTransformer{<:ZScoreTransform}, GaussianProcessInputTransformer{<:UnitRangeTransform}} +) = to_gp_format( + StatsBase.transform( + transformer.transform, + dataframe_to_array(data, transformer.input) + ) +) + +function transform( + data::DataFrame, + transformer::GaussianProcessInputTransformer{<:StandardNormalTransform} + ) + data_copy = copy(data) + to_standard_normal_space!(transformer.input, data_copy) + return to_gp_format(dataframe_to_array(data_copy, names(transformer.input))) +end + + +# ---------------- Output transformation ---------------- +# # Developer Note +# Gaussian process regression target outputs are always extracted from DataFrame and transformed to Vector +# Inverse transforms return Vectors that later get inserted in the provided DataFrame that is used in evaluate! method for GaussianProcess +struct GaussianProcessOutputTransformer{T} + transform::T + output::Symbol +end + +# Fitting +fit_output_transform( + ::DataFrame, + output::Symbol, + ::IdentityTransformChoice +) = GaussianProcessOutputTransformer(NoTransform(), output) + +function fit_output_transform( + data::DataFrame, + output::Symbol, + ::ZScoreTransformChoice + ) + transform = fit( + StatsBase.ZScoreTransform, + dataframe_to_array(data, output); + dims = 1 + ) + return GaussianProcessOutputTransformer(transform, output) +end + +function fit_output_transform( + data::DataFrame, + output::Symbol, + ::UnitRangeTransformChoice + ) + transform = fit( + StatsBase.UnitRangeTransform, + dataframe_to_array(data, output); + dims = 1 + ) + return GaussianProcessOutputTransformer(transform, output) +end + +fit_output_transform( + ::DataFrame, + ::Symbol, + ::StandardNormalTransformChoice +) = throw(ArgumentError("Standard normal transform for outputs is not possible")) + +# Transforms +transform( + data::DataFrame, + transformer::GaussianProcessOutputTransformer{NoTransform} +) = to_gp_format(dataframe_to_array(data, transformer.output)) + +transform( + data::DataFrame, + transformer::Union{GaussianProcessOutputTransformer{<:ZScoreTransform}, GaussianProcessOutputTransformer{<:UnitRangeTransform}} +) = to_gp_format( + StatsBase.transform( + transformer.transform, + dataframe_to_array(data, transformer.output) + ) +) + +# Inverse transforms +# # Developer Note +# Gaussian process regression requires two distinct inverse transformations for the output: +# one for the mean predictions (this same transformation can also be applied to function samples) and one for the variance predictions. + +# Consider a z-score transformation of output ``y``: +# ```math +# \tilde{y} = \frac{y - μ}{σ}. +# ``` +# To recover the mean of the untransformed output, we can simply apply the inverse transformation: +# ```math +# E[y] = E[σ\tilde{y} + μ] = σE[\tilde{y}] + μ. +# ``` +# Analogously, sampled functions ``\tilde{y}_s`` from the Gaussian process regression model can be transformed back: +# ```math +# y_s = σ\tilde{y}_s + μ. +# ``` +# The variance, however, is untransformed as follows: +# ```math +# Var[y] = E[(σ\tilde{y} + μ - E[σ\tilde{y} + μ])^2] = E[(σ^2(\tilde{y} - E[\tilde{y}])^2] = σ^2 Var[\tilde{y}] +# ``` +inverse_transform( + data::AbstractArray, + ::GaussianProcessOutputTransformer{NoTransform} +) = data + +inverse_transform( + data::AbstractArray, + transformer::Union{GaussianProcessOutputTransformer{<:ZScoreTransform}, GaussianProcessOutputTransformer{<:UnitRangeTransform}} +) = StatsBase.reconstruct(transformer.transform, data) + +variance_inverse_transform( + data::AbstractArray, + ::GaussianProcessOutputTransformer{NoTransform} +) = data + +variance_inverse_transform( + data::AbstractArray, + transformer::Union{GaussianProcessOutputTransformer{<:ZScoreTransform}, GaussianProcessOutputTransformer{<:UnitRangeTransform}} +) = only(transformer.transform.scale)^2 * data diff --git a/test/Project.toml b/test/Project.toml index b93de7fe6..dfa9cc34e 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,11 +1,15 @@ [deps] +AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918" Copulas = "ae264745-0b69-425e-9d9d-cf662c5eec93" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +ParameterHandling = "2412ca09-6db7-441c-8e3a-88d5709968c5" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/models/gp/adaptivegaussianprocess.jl b/test/models/gp/adaptivegaussianprocess.jl new file mode 100644 index 000000000..d228beb87 --- /dev/null +++ b/test/models/gp/adaptivegaussianprocess.jl @@ -0,0 +1,123 @@ +@testset "Adaptive Gaussian Process" begin + + # Input + input = RandomVariable(Uniform(-2, 12), :x1) + n_design_points = 8 + n_added_points = 3 + design = LatinHypercubeSampling(n_design_points) + model = Model(df -> df.x1 .^ 2 .* sin.(df.x1), :y) + + acquisition_function = MaximumVariance() + # small candidate set keeps the tests fast + candidate_sampling = MonteCarlo(200) + + prior = GP(ConstMean(0.0), SqExponentialKernel()) + + # Initial design used by the DataFrame-based constructors + data = sample(input, design) + evaluate!(model, data) + + @testset "GP prior + UQInput" begin + + gp = AdaptiveGaussianProcess( + prior, input, model, :y, acquisition_function, n_added_points, + n_design_points, design; + candidate_sampling = candidate_sampling, + ) + + @test gp isa GaussianProcess + @test size(gp.training_data, 1) == n_design_points + n_added_points + @test gp.output == :y + end + + @testset "Default prior + UQInput" begin + + gp_default = AdaptiveGaussianProcess( + input, model, :y, acquisition_function, n_added_points, + n_design_points, design; + candidate_sampling = candidate_sampling, + ) + + @test gp_default isa GaussianProcess + @test size(gp_default.training_data, 1) == n_design_points + n_added_points + end + + @testset "Pre-fit GaussianProcess" begin + gp_model = GaussianProcess( + input, model, :y; + experimental_design = design, + mean_fct = ZeroMean(), + kernel = SqExponentialKernel(), + ) + + gp = AdaptiveGaussianProcess( + gp_model, input, model, acquisition_function, n_added_points; + candidate_sampling = candidate_sampling, + ) + + @test gp isa GaussianProcess + @test size(gp.training_data, 1) == n_design_points + n_added_points + + # n_added_points = 0 should return the training data unchanged + gp_unchanged = AdaptiveGaussianProcess( + gp_model, input, model, acquisition_function, 0; + candidate_sampling = candidate_sampling, + ) + @test size(gp_unchanged.training_data, 1) == n_design_points + n_added_points + end + + @testset "GP prior + DataFrame" begin + + gp = AdaptiveGaussianProcess( + prior, copy(data), input, model, :y, acquisition_function, n_added_points; + candidate_sampling = candidate_sampling, + ) + + @test gp isa GaussianProcess + @test size(gp.training_data, 1) == n_design_points + n_added_points + @test gp.output == :y + end + + @testset "Default prior + DataFrame" begin + + gp_default = AdaptiveGaussianProcess( + copy(data), input, model, :y, acquisition_function, n_added_points; + candidate_sampling = candidate_sampling, + ) + + gp_explicit = AdaptiveGaussianProcess( + prior, copy(data), input, model, :y, acquisition_function, n_added_points; + candidate_sampling = candidate_sampling, + ) + + @test gp_default isa GaussianProcess + @test size(gp_default.training_data, 1) == n_design_points + n_added_points + end + + @testset "Not learn hyperparameters" begin + gp_model = GaussianProcess( + input, model, :y; + experimental_design = design, + mean_fct = ConstMean(0.0), + kernel = MaternKernel() + ) + + new_data = sample(input) + evaluate!(model, new_data) + + new_gp_model = UncertaintyQuantification._refit_gp( + deepcopy(gp_model), new_data, MaximumLikelihoodEstimation(), 1.0e-9, false, false + ) + + trend_initial = gp_model.posterior.prior.mean.c + trend_adaptive = new_gp_model.posterior.prior.mean.c + + kernel_initial = gp_model.posterior.prior.kernel.ν + kernel_adaptive = new_gp_model.posterior.prior.kernel.ν + + @test trend_initial == trend_adaptive + @test kernel_initial == kernel_adaptive + + end + +end diff --git a/test/models/gp/gaussianprocess.jl b/test/models/gp/gaussianprocess.jl new file mode 100644 index 000000000..d50ddf1a3 --- /dev/null +++ b/test/models/gp/gaussianprocess.jl @@ -0,0 +1,352 @@ +function create_test_data(n_samples::Int, lower::Real, upper::Real, dim::Int) + data = lower .+ (upper - lower) .* rand(n_samples, dim) + df = DataFrame() + for i in 1:dim + name = Symbol("x$i") + df[!, name] = data[:, i] + end + return df +end + +@testset "Gaussian Process" begin + # Input samples + n_input_samples = 10 + design = LatinHypercubeSampling(n_input_samples) + # Lower and upper bound for inputs + lower = 0 + upper = 5 + + # Use same base gp for every test + σ² = 1.0e-5 + mean_fct = ConstMean(0.0) + kernel = SqExponentialKernel() + + # Possible transforms + input_transform_choices = [ + IdentityTransformChoice, StandardNormalTransformChoice, + UnitRangeTransformChoice, ZScoreTransformChoice, + ] + output_transform_choices = [ + IdentityTransformChoice, UnitRangeTransformChoice, ZScoreTransformChoice, + ] + + @testset "GP Optimization" begin + @testset "Dataframe" begin + x = collect(range(lower, stop = upper, length = n_input_samples)) + y = sin.(x) + data = DataFrame(:x1 => x, :y => y) + gp = GaussianProcess( + data, :y; + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = 0.0, + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + #check if σ² was left unchanged in the optimization + @test gp.σ² == 0.0 + + gp = GaussianProcess( + data, :y; + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = 0.0, + mean_fct = mean_fct, + kernel = kernel, + learn_noise = true + ) + @test gp.σ² > 0.0 + + gp = GaussianProcess( + data, :y; + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + @test gp.σ² == σ² + + gp = GaussianProcess( + data, :y; + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = true + ) + @test gp.σ² != σ² + + df_mean = create_test_data(n_input_samples, lower, upper, 1) + df_var = create_test_data(n_input_samples, lower, upper, 1) + df_mean_var = create_test_data(n_input_samples, lower, upper, 1) + df_samples = create_test_data(n_input_samples, lower, upper, 1) + + evaluate!(gp, df_mean; mode = :mean) + evaluate!(gp, df_var; mode = :var) + evaluate!(gp, df_mean_var; mode = :mean_and_var) + evaluate!(gp, df_samples; mode = :sample, n_samples = 1) + @test :y_mean in propertynames(df_mean) + @test !(:y_var in propertynames(df_mean)) + @test :y_var in propertynames(df_var) + @test !(:y_mean in propertynames(df_var)) + @test :y_mean in propertynames(df_mean_var) + @test :y_var in propertynames(df_mean_var) + @test :y_sample_1 in propertynames(df_samples) + @test_throws ArgumentError evaluate!(gp, df_mean; mode = :error) + + @test_throws DomainError gp = GaussianProcess( + data, :y; + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = -1.0, + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + end + @testset "UQInput" begin + xrv = [Parameter(1.5, :p), RandomVariable(Uniform(lower, upper), :x1)] + xrv_single = RandomVariable(Uniform(lower, upper), :x1) + model = Model( + df -> df.p .* sin.(df.x1), :y + ) + model_single = Model( + df -> 1.5 .* sin.(df.x1), :y + ) + + gp = GaussianProcess( + xrv_single, model_single, :y; + experimental_design = design, + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = 0.0, + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + + gp = GaussianProcess( + xrv, model, :y; + experimental_design = design, + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = 0.0, + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + @test gp.σ² == 0.0 + + gp = GaussianProcess( + xrv, model, :y; + experimental_design = design, + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = 0.0, + mean_fct = mean_fct, + kernel = kernel, + learn_noise = true + ) + @test gp.σ² > 0.0 + + gp = GaussianProcess( + xrv, model, :y; + experimental_design = design, + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + @test gp.σ² == σ² + + gp = GaussianProcess( + xrv, model, :y; + experimental_design = design, + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = true + ) + @test gp.σ² != σ² + + df_mean = sample(xrv, n_input_samples) + df_var = sample(xrv, n_input_samples) + df_mean_var = sample(xrv, n_input_samples) + evaluate!(gp, df_mean; mode = :mean) + evaluate!(gp, df_var; mode = :var) + evaluate!(gp, df_mean_var; mode = :mean_and_var) + @test :y_mean in propertynames(df_mean) + @test !(:y_var in propertynames(df_mean)) + @test :y_var in propertynames(df_var) + @test !(:y_mean in propertynames(df_var)) + @test :y_mean in propertynames(df_mean_var) + @test :y_var in propertynames(df_mean_var) + + @test_throws DomainError gp = GaussianProcess( + xrv, model, :y; + experimental_design = design, + input_transform = IdentityTransformChoice(), + output_transform = IdentityTransformChoice(), + σ² = -1.0, + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + end + end + @testset "1D Input" begin + @testset "Dataframe" begin + x = collect(range(lower, stop = upper, length = n_input_samples)) + y = sin.(x) + data = DataFrame(:x1 => x, :y => y) + for input_transform in input_transform_choices, output_transform in output_transform_choices + @testset "$input_transform → $output_transform" begin + if input_transform == StandardNormalTransformChoice + @test_throws ArgumentError GaussianProcess( + data, :y; + input_transform = input_transform(), + output_transform = output_transform(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel + ) + else + gp = GaussianProcess( + data, :y; + input_transform = input_transform(), + output_transform = output_transform(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + #check if σ² was left unchanged in the optimization + @test gp.σ² == σ² + df = create_test_data(n_input_samples, lower, upper, 1) + evaluate!(gp, df; mode = :mean_and_var) + @test :y_mean in propertynames(df) + @test :y_var in propertynames(df) + end + end + end + end + @testset "UQInput" begin + xrv = [Parameter(1.5, :p), RandomVariable(Uniform(lower, upper), :x1)] + model = Model( + df -> df.p .* sin.(df.x1), :y + ) + for input_transform in input_transform_choices, output_transform in output_transform_choices + @testset "$input_transform → $output_transform" begin + if input_transform == StandardNormalTransformChoice + @test_throws ArgumentError GaussianProcess( + xrv, model, :y; + experimental_design = design, + input_transform = input_transform(), + output_transform = output_transform(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + else + gp = GaussianProcess( + xrv, model, :y; + experimental_design = design, + input_transform = input_transform(), + output_transform = output_transform(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + df = sample(xrv, n_input_samples) + evaluate!(gp, df; mode = :mean_and_var) + @test :y_mean in propertynames(df) + @test :y_var in propertynames(df) + end + end + end + end + end + @testset "2D Input" begin + @testset "Dataframe" begin + x = [collect(range(lower, stop = upper, length = n_input_samples)) collect(range(lower, stop = upper, length = n_input_samples))] + y = sin.(x[:, 1]) + cos.(x[:, 2]) + data = DataFrame(:x1 => x[:, 1], :x2 => x[:, 2], :y => y) + for input_transform in input_transform_choices, output_transform in output_transform_choices + @testset "$input_transform → $output_transform" begin + if input_transform == StandardNormalTransformChoice + @test_throws ArgumentError GaussianProcess( + data, :y; + input_transform = input_transform(), + output_transform = output_transform(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + else + gp = GaussianProcess( + data, :y; + input_transform = input_transform(), + output_transform = output_transform(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + df = create_test_data(n_input_samples, lower, upper, 2) + evaluate!(gp, df; mode = :mean_and_var) + @test :y_mean in propertynames(df) + @test :y_var in propertynames(df) + end + end + end + end + @testset "UQInput" begin + xrv = [Parameter(1.5, :p), RandomVariable(Uniform(0, 5), :x1), RandomVariable(Uniform(0, 5), :x2)] + model = Model( + df -> df.p .* sin.(df.x1) + df.p .* cos.(df.x2), :y + ) + for input_transform in input_transform_choices, output_transform in output_transform_choices + @testset "$input_transform → $output_transform" begin + if input_transform == StandardNormalTransformChoice + @test_throws ArgumentError GaussianProcess( + xrv, model, :y; + experimental_design = design, + input_transform = input_transform(), + output_transform = output_transform(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + else + gp = GaussianProcess( + xrv, model, :y; + experimental_design = design, + input_transform = input_transform(), + output_transform = output_transform(), + σ² = σ², + mean_fct = mean_fct, + kernel = kernel, + learn_noise = false + ) + df = sample(xrv, n_input_samples) + evaluate!(gp, df; mode = :mean_and_var) + @test :y_mean in propertynames(df) + @test :y_var in propertynames(df) + end + end + end + end + end +end diff --git a/test/models/gp/gp_acquisitionfunction.jl b/test/models/gp/gp_acquisitionfunction.jl new file mode 100644 index 000000000..c166d414d --- /dev/null +++ b/test/models/gp/gp_acquisitionfunction.jl @@ -0,0 +1,141 @@ +@testset "GP Acquisition Functions" begin + # Small 1D training set: y = x^2 * sin(x) + input = RandomVariable(Uniform(-2, 12), :x1) + n_input_samples = 10 + design = LatinHypercubeSampling(n_input_samples) + model = Model(df -> df.x1 .^ 2 .* sin.(df.x1), :y) + + data = sample(input, design) + evaluate!(model, data) + + gp = GaussianProcess( + data, :y; + mean_fct = ConstMean(0.0), + kernel = SqExponentialKernel(), + learn_hyperparameters = true, + ) + + candidates = DataFrame(x1 = collect(range(-2, 12, 25))) + + function mean_and_std(gp, candidates) + df = copy(candidates) + evaluate!(gp, df; mode = :mean_and_var) + return df[:, :y_mean], sqrt.(df[:, :y_var]) + end + + @testset "MaximumVariance" begin + μ, σ = mean_and_std(gp, candidates) + expected = candidates[[argmax(σ .^ 2)], [:x1]] + + result = UncertaintyQuantification._find_next_point(gp, candidates, MaximumVariance()) + @test result == expected + end + + @testset "ExpectedImprovement" begin + ei = ExpectedImprovement(; ξ = 0.0) + μ, σ = mean_and_std(gp, candidates) + f_best = minimum(data.y) + + z = (f_best .- μ) ./ σ + ei_values = (f_best .- μ) .* cdf.(Normal(), z) .+ σ .* pdf.(Normal(), z) + expected = candidates[[argmax(ei_values)], [:x1]] + + result = UncertaintyQuantification._find_next_point(gp, candidates, ei) + @test result == expected + end + + @testset "ProbabilityOfImprovement" begin + poi = ProbabilityOfImprovement(; ξ = 0.0) + μ, σ = mean_and_std(gp, candidates) + f_best = minimum(data.y) + + z = (f_best .- μ) ./ σ + pi_values = cdf.(Normal(), z) + expected = candidates[[argmax(pi_values)], [:x1]] + + result = UncertaintyQuantification._find_next_point(gp, candidates, poi) + @test result == expected + end + + @testset "UpperConfidenceBound" begin + ucb = UpperConfidenceBound(; κ = 2.0) + μ, σ = mean_and_std(gp, candidates) + + lcb = μ .- ucb.κ .* σ + expected = candidates[[argmin(lcb)], [:x1]] + + result = UncertaintyQuantification._find_next_point(gp, candidates, ucb) + @test result == expected + + # larger κ should not decrease the selected lower confidence bound value + ucb_explore = UpperConfidenceBound(; κ = 10.0) + result_explore = UncertaintyQuantification._find_next_point(gp, candidates, ucb_explore) + @test result_explore isa DataFrame + end + + @testset "DeviationNumber" begin + dn = DeviationNumber(; threshold = 0.0, stopping = 2.0) + μ, σ = mean_and_std(gp, candidates) + + u = abs.(μ .- dn.threshold) ./ σ + expected = candidates[[argmin(u)], [:x1]] + + result = UncertaintyQuantification._find_next_point(gp, candidates, dn) + @test result == expected + + # stopping criterion is true once min(U) exceeds the threshold + _, converged = UncertaintyQuantification._find_next_point_stopping(gp, candidates, dn) + @test converged == (minimum(u) >= dn.stopping) + end + + @testset "ExpectedFeasibility" begin + eff = ExpectedFeasibility(; threshold = 0.0, epsilon_factor = 2.0, stopping = 0.001) + μ, σ = mean_and_std(gp, candidates) + ε = eff.epsilon_factor .* σ + + z = (eff.threshold .- μ) ./ σ + z_lower = (eff.threshold .- ε .- μ) ./ σ + z_upper = (eff.threshold .+ ε .- μ) ./ σ + + eff_values = + (μ .- eff.threshold) .* + (2 .* cdf.(Normal(), z) .- cdf.(Normal(), z_lower) .- cdf.(Normal(), z_upper)) .- + σ .* (2 .* pdf.(Normal(), z) .- pdf.(Normal(), z_lower) .- pdf.(Normal(), z_upper)) .+ + ε .* (cdf.(Normal(), z_upper) .- cdf.(Normal(), z_lower)) + + expected = candidates[[argmax(eff_values)], [:x1]] + + result = UncertaintyQuantification._find_next_point(gp, candidates, eff) + @test result == expected + + _, converged = UncertaintyQuantification._find_next_point_stopping(gp, candidates, eff) + @test converged == (maximum(eff_values) <= eff.stopping) + end + + @testset "MaximinDistance" begin + # doesn't depend on the GP posterior: purely geometric + X = Matrix(data[:, [:x1]]) + Xc = Matrix(candidates[:, [:x1]]) + distances = [minimum(norm(Xc[i, :] - X[j, :]) for j in axes(X, 1)) for i in axes(Xc, 1)] + expected = candidates[[argmax(distances)], [:x1]] + + result = UncertaintyQuantification._find_next_point(gp, candidates, MaximinDistance()) + @test result == expected + end + + @testset "ExpectedImprovementForGlobalFit" begin + μ, σ = mean_and_std(gp, candidates) + σ² = σ .^ 2 + + X = Matrix(data[:, [:x1]]) + Xc = Matrix(candidates[:, [:x1]]) + y = data.y + + nearest = [argmin([norm(Xc[i, :] - X[j, :]) for j in axes(X, 1)]) for i in axes(Xc, 1)] + eigf = abs2.(μ .- y[nearest]) .+ σ² + expected = candidates[[argmax(eigf)], [:x1]] + + result = UncertaintyQuantification._find_next_point(gp, candidates, ExpectedImprovementForGlobalFit()) + @test result == expected + end +end diff --git a/test/models/gp/hyperparametertuning.jl b/test/models/gp/hyperparametertuning.jl new file mode 100644 index 000000000..a2b068bee --- /dev/null +++ b/test/models/gp/hyperparametertuning.jl @@ -0,0 +1,39 @@ +@testset "Hyperparameter tuning" begin + + @testset "1D Input" begin + x = collect(range(0, stop = 10, length = 5)) + fct(x) = sin.(x) .+ x + y = fct(x) + data = DataFrame(:x => x, :y => y) + + σ² = 1.0e-9 + kernel = SqExponentialKernel() ∘ ScaleTransform(10.0) + prior_gp = GP(ConstMean(0.0), kernel) + gp_opt = GaussianProcess(prior_gp, data, :y; σ² = σ²) + gp_nonopt = GaussianProcess(prior_gp, data, :y; σ² = σ², learn_hyperparameters = false) + x_test = collect(range(0, stop = 10, length = 50)) + y_test = fct(x_test) + likelihood_no_opt = logpdf(gp_nonopt.posterior(x_test), y_test) + likelihood_opt = logpdf(gp_opt.posterior(x_test), y_test) + + @test likelihood_opt > likelihood_no_opt + end + + @testset "2D Input" begin + x = [collect(range(0, stop = 5, length = 10)) collect(range(0, stop = 5, length = 10))] + y = sin.(x[:, 1]) + cos.(x[:, 2]) + data = DataFrame(:x1 => x[:, 1], :x2 => x[:, 2], :y => y) + + σ² = 1.0e-9 + kernel = Matern52Kernel() ∘ ARDTransform([5.0, 5.0]) + prior_gp = GP(ConstMean(0.0), kernel) + gp_opt = GaussianProcess(prior_gp, data, :y; σ² = σ²) + gp_nonopt = GaussianProcess(prior_gp, data, :y; σ² = σ², learn_hyperparameters = false) + + x_test = [collect(range(0, stop = 5, length = 50)) collect(range(0, stop = 5, length = 50))] + y_test = sin.(x_test[:, 1]) + cos.(x_test[:, 2]) + likelihood_no_opt = logpdf(gp_nonopt.posterior(RowVecs(x_test)), y_test) + likelihood_opt = logpdf(gp_opt.posterior(RowVecs(x_test)), y_test) + @test likelihood_opt > likelihood_no_opt + end +end diff --git a/test/models/gp/parameterization.jl b/test/models/gp/parameterization.jl new file mode 100644 index 000000000..2b808a4df --- /dev/null +++ b/test/models/gp/parameterization.jl @@ -0,0 +1,183 @@ +function is_of_type(exporting_module::Module, name::Symbol, type::DataType) + obj = getfield(exporting_module, name) + if obj isa DataType + return obj <: type + elseif obj isa UnionAll + return obj.body <: type + else + return false + end +end + +function get_exported_types(exporting_module::Module, type::DataType) + exported_names = names(exporting_module; all = false) + type_symbols = filter(n -> is_of_type(exporting_module, n, type), exported_names) + types = map(sym -> getfield(exporting_module, sym), type_symbols) + return filter(t -> !isabstracttype(t), types) +end + +function check_extract_parameters(type::Type) + return hasmethod(UncertaintyQuantification.extract_parameters, Tuple{type}) +end +function check_apply_parameters(type::Type) + return hasmethod(UncertaintyQuantification.apply_parameters, Tuple{type, Any}) +end +function check_implementation(type::Type) + return check_extract_parameters(type) && check_apply_parameters(type) +end + +@testset "Parameterization" begin + @testset "Mean functions" begin + @testset "ZeroMean" begin + m = ZeroMean() + @test isnothing(UncertaintyQuantification.extract_parameters(m)) + @test UncertaintyQuantification.apply_parameters(m, nothing) === m + end + @testset "ConstMean" begin + m = ConstMean(2.5) + θ = UncertaintyQuantification.extract_parameters(m) + @test θ ≈ 2.5 + @test UncertaintyQuantification.apply_parameters(m, θ).c ≈ 2.5 + end + @testset "CustomMean" begin + m = CustomMean(x -> sum(x)) + @test isnothing(UncertaintyQuantification.extract_parameters(m)) + @test UncertaintyQuantification.apply_parameters(m, nothing) === m + end + @testset "Unimplemented Means" begin + # Check if any means exported from KernelFunctions.jl are not handled by parameterization + meanfunctions = get_exported_types(AbstractGPs, AbstractGPs.MeanFunction) + unimplemented_meanfunctions = filter(!check_implementation, meanfunctions) + if !isempty(unimplemented_meanfunctions) + @error "Mean parameter handling not implemented for:\n " * + join(string.(unimplemented_meanfunctions), "\n ") + end + @test isempty(unimplemented_meanfunctions) + end + end + + @testset "Kernel functions" begin + # @testset "Kernels without parameters" begin + # TODO: Test kernels without parameters + # end + + @testset "Kernels with parameters" begin + @testset "ConstantKernel" begin + k = ConstantKernel(; c = 3.0) + θ = UncertaintyQuantification.extract_parameters(k) + k2 = UncertaintyQuantification.apply_parameters( + k, ParameterHandling.value(θ) + ) + @test only(k2.c) ≈ 3.0 + end + @testset "ScaleTransform" begin + t = ScaleTransform(2.0) + θ = UncertaintyQuantification.extract_parameters(t) + t2 = UncertaintyQuantification.apply_parameters( + t, ParameterHandling.value(θ) + ) + @test only(t2.s) ≈ 2.0 + end + @testset "ARDTransform" begin + t = ARDTransform([1.0, 2.0, 3.0]) + θ = UncertaintyQuantification.extract_parameters(t) + t2 = UncertaintyQuantification.apply_parameters( + t, ParameterHandling.value(θ) + ) + @test t2.v ≈ [1.0, 2.0, 3.0] + end + @testset "PeriodicKernel" begin + k = PeriodicKernel(; r = [1.5]) + θ = UncertaintyQuantification.extract_parameters(k) + k2 = UncertaintyQuantification.apply_parameters( + k, ParameterHandling.value(θ) + ) + @test only(k2.r) ≈ 1.5 + end + @testset "ScaledKernel" begin + k = 2.0 * SqExponentialKernel() # ScaledKernel + θ = UncertaintyQuantification.extract_parameters(k) + k2 = UncertaintyQuantification.apply_parameters( + k, ParameterHandling.value(θ) + ) + @test only(k2.σ²) ≈ 2.0 + end + @testset "RationalQuadraticKernel" begin + k = RationalQuadraticKernel(; α = 1.5) + θ = UncertaintyQuantification.extract_parameters(k) + k2 = UncertaintyQuantification.apply_parameters( + k, ParameterHandling.value(θ) + ) + @test only(k2.α) ≈ 1.5 + end + end + + @testset "Composite kernels" begin + @testset "KernelSum" begin + k = SqExponentialKernel() + ConstantKernel(; c = 2.0) + θ = UncertaintyQuantification.extract_parameters(k) + k2 = UncertaintyQuantification.apply_parameters( + k, ParameterHandling.value(θ) + ) + # first component has no params, second has c=2.0 + @test isnothing(θ[1]) + @test only(k2.kernels[2].c) ≈ 2.0 + end + @testset "KernelProduct" begin + k = SqExponentialKernel() * ConstantKernel(; c = 3.0) + θ = UncertaintyQuantification.extract_parameters(k) + k2 = UncertaintyQuantification.apply_parameters( + k, ParameterHandling.value(θ) + ) + @test only(k2.kernels[2].c) ≈ 3.0 + end + @testset "TransformedKernel" begin + k = SqExponentialKernel() ∘ ScaleTransform(2.0) + θ = UncertaintyQuantification.extract_parameters(k) + k2 = UncertaintyQuantification.apply_parameters( + k, ParameterHandling.value(θ) + ) + @test only(k2.transform.s) ≈ 2.0 + end + @testset "Nested composite: ScaledKernel with TransformedKernel" begin + k = 4.0 * (SqExponentialKernel() ∘ ScaleTransform(2.0)) + θ = UncertaintyQuantification.extract_parameters(k) + k2 = UncertaintyQuantification.apply_parameters( + k, ParameterHandling.value(θ) + ) + @test only(k2.σ²) ≈ 4.0 + @test only(k2.kernel.transform.s) ≈ 2.0 + end + end + + @testset "Unimplemented Kernels" begin + # Check if any kernels and transforms exported from KernelFunctions.jl are not handled by parameterization + transforms = get_exported_types(KernelFunctions, Transform) + unimplemented_transforms = filter(!check_implementation, transforms) + kernels = get_exported_types(KernelFunctions, Kernel) + unimplemented_kernels = filter(!check_implementation, kernels) + + if !isempty(unimplemented_transforms) + @error "Transform parameter handling not implemented for:\n " * + join(string.(unimplemented_transforms), "\n ") + end + @test isempty(unimplemented_transforms) + + if !isempty(unimplemented_kernels) + @error "Kernel parameter handling not implemented for:\n " * + join(string.(unimplemented_kernels), "\n ") + end + @test isempty(unimplemented_kernels) + end + end + + @testset "GP" begin + @testset "GP round-trip" begin + gp = GP(ConstMean(1.0), SqExponentialKernel()) + θ = UncertaintyQuantification.extract_parameters(gp) + gp2 = UncertaintyQuantification.apply_parameters(gp, ParameterHandling.value(θ)) + @test only(gp2.mean.c) ≈ 1.0 + @test gp2.kernel isa SqExponentialKernel + end + end +end diff --git a/test/models/gp/standardization.jl b/test/models/gp/standardization.jl new file mode 100644 index 000000000..9e7c7aaa1 --- /dev/null +++ b/test/models/gp/standardization.jl @@ -0,0 +1,122 @@ +@testset "Standardization" begin + + @testset "to_gp_format" begin + v = [1.0, 2.0, 3.0] + @test UncertaintyQuantification.to_gp_format(v) isa Vector + M = [1.0 2.0; 3.0 4.0; 5.0 6.0] + @test UncertaintyQuantification.to_gp_format(M) isa RowVecs + end + + @testset "Input transformer" begin + @testset "IdentityTransformChoice" begin + data = DataFrame(x = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_input_transform(data, :x, IdentityTransformChoice()) + x = UncertaintyQuantification.transform(data, transformer) + @test x ≈ data.x + end + @testset "ZScoreTransformChoice" begin + data = DataFrame(x = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_input_transform(data, :x, ZScoreTransformChoice()) + x = UncertaintyQuantification.transform(data, transformer) + @test mean(x) ≈ 0.0 atol = 1.0e-10 + @test std(x) ≈ 1.0 atol = 1.0e-10 + end + @testset "UnitRangeTransformChoice" begin + data = DataFrame(x = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_input_transform(data, :x, UnitRangeTransformChoice()) + x = UncertaintyQuantification.transform(data, transformer) + @test minimum(x) ≈ 0.0 atol = 1.0e-10 + @test maximum(x) ≈ 1.0 atol = 1.0e-10 + end + @testset "StandardNormalTransformChoice with Symbol throws" begin + data = DataFrame(x = [1.0, 2.0, 3.0]) + @test_throws ArgumentError UncertaintyQuantification.fit_input_transform(data, :x, StandardNormalTransformChoice()) + end + @testset "Multivariate input gives RowVecs" begin + data = DataFrame(x1 = [1.0, 2.0, 3.0], x2 = [4.0, 5.0, 6.0]) + transformer = UncertaintyQuantification.fit_input_transform(data, [:x1, :x2], IdentityTransformChoice()) + x = UncertaintyQuantification.transform(data, transformer) + @test x isa RowVecs + end + @testset "Single input column gives Vector" begin + data = DataFrame(x = [1.0, 2.0, 3.0]) + transformer = UncertaintyQuantification.fit_input_transform(data, :x, IdentityTransformChoice()) + x = UncertaintyQuantification.transform(data, transformer) + @test x isa Vector + end + end + + @testset "Output transformer" begin + @testset "IdentityTransformChoice" begin + data = DataFrame(y = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_output_transform(data, :y, IdentityTransformChoice()) + y = UncertaintyQuantification.transform(data, transformer) + @test y ≈ data.y + end + @testset "ZScoreTransformChoice" begin + data = DataFrame(y = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_output_transform(data, :y, ZScoreTransformChoice()) + y = UncertaintyQuantification.transform(data, transformer) + @test mean(y) ≈ 0.0 atol = 1.0e-10 + @test std(y) ≈ 1.0 atol = 1.0e-10 + end + @testset "UnitRangeTransformChoice" begin + data = DataFrame(y = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_output_transform(data, :y, UnitRangeTransformChoice()) + y = UncertaintyQuantification.transform(data, transformer) + @test minimum(y) ≈ 0.0 atol = 1.0e-10 + @test maximum(y) ≈ 1.0 atol = 1.0e-10 + end + @testset "StandardNormalTransformChoice throws" begin + data = DataFrame(y = [1.0, 2.0, 3.0]) + @test_throws ArgumentError UncertaintyQuantification.fit_output_transform(data, :y, StandardNormalTransformChoice()) + end + end + + @testset "Inverse transforms" begin + @testset "NoTransform is identity" begin + data = DataFrame(y = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_output_transform(data, :y, IdentityTransformChoice()) + y = UncertaintyQuantification.transform(data, transformer) + @test UncertaintyQuantification.inverse_transform(y, transformer) ≈ data.y + @test UncertaintyQuantification.variance_inverse_transform(y, transformer) ≈ y + end + @testset "ZScoreTransform round-trip" begin + data = DataFrame(y = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_output_transform(data, :y, ZScoreTransformChoice()) + y_transformed = UncertaintyQuantification.transform(data, transformer) + @test UncertaintyQuantification.inverse_transform(y_transformed, transformer) ≈ data.y atol = 1.0e-10 + end + @testset "UnitRangeTransform round-trip" begin + data = DataFrame(y = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_output_transform(data, :y, UnitRangeTransformChoice()) + y_transformed = UncertaintyQuantification.transform(data, transformer) + @test UncertaintyQuantification.inverse_transform(y_transformed, transformer) ≈ data.y atol = 1.0e-10 + end + @testset "Variance inverse transform - ZScore" begin + data = DataFrame(y = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_output_transform(data, :y, ZScoreTransformChoice()) + σ² = only(transformer.transform.scale)^2 + var_transformed = ones(5) + # Var[y] = σ² * Var[ỹ] + @test UncertaintyQuantification.variance_inverse_transform(var_transformed, transformer) ≈ σ² * var_transformed atol = 1.0e-10 + end + @testset "Variance inverse transform - UnitRange" begin + data = DataFrame(y = [1.0, 2.0, 3.0, 4.0, 5.0]) + transformer = UncertaintyQuantification.fit_output_transform(data, :y, UnitRangeTransformChoice()) + σ² = only(transformer.transform.scale)^2 + var_transformed = ones(5) + @test UncertaintyQuantification.variance_inverse_transform(var_transformed, transformer) ≈ σ² * var_transformed atol = 1.0e-10 + end + @testset "Variance inverse is not the same as mean inverse" begin + # Regression guard: make sure variance and mean inverse transforms are not accidentally swapped + data = DataFrame(y = [10.0, 20.0, 30.0, 40.0, 50.0]) + transformer = UncertaintyQuantification.fit_output_transform(data, :y, ZScoreTransformChoice()) + var_transformed = ones(5) + mean_inv = UncertaintyQuantification.inverse_transform(var_transformed, transformer) + var_inv = UncertaintyQuantification.variance_inverse_transform(var_transformed, transformer) + @test !(mean_inv ≈ var_inv) + end + end + +end diff --git a/test/runtests.jl b/test/runtests.jl index aba99e9ef..93b050396 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,12 +4,15 @@ using Distributed using HCubature using HypothesisTests using InteractiveUtils +using LinearAlgebra: norm +using ParameterHandling using Random using StatsBase: fit, Histogram, corkendall using Test using Plots using TestItemRunner using UncertaintyQuantification +import UncertaintyQuantification: sample @testsnippet TestSetup begin using Copulas @@ -46,15 +49,21 @@ end end end + include("models/model.jl") +include("models/gp/gaussianprocess.jl") +include("models/gp/hyperparametertuning.jl") +include("models/gp/parameterization.jl") +include("models/gp/standardization.jl") +include("models/gp/gp_acquisitionfunction.jl") +include("models/gp/adaptivegaussianprocess.jl") include("modelupdating/bayesianTM.jl") include("inputs/jointdistribution.jl") include("inputs/imprecise/p-box.jl") -@run_package_tests - - include("plotting/plotting.jl") +@run_package_tests + if Sys.islinux() HPC = false HPC_account = "HPC_account_1"