From f0d089671b5e08183c29cb10322d2090f45fefc3 Mon Sep 17 00:00:00 2001 From: odunbar Date: Mon, 29 Jun 2020 17:08:50 +0100 Subject: [PATCH 01/18] GPEmulator.jl added SVD functionality --- src/GPEmulator.jl | 238 ++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 232 insertions(+), 6 deletions(-) diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 834e74e0f..12fa5de03 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -12,7 +12,7 @@ using PyCall @sk_import gaussian_process : GaussianProcessRegressor @sk_import gaussian_process.kernels : (RBF, WhiteKernel, ConstantKernel) -export GPObj +export GPObj, NormalizedSVDGPObj export predict export GPJL, SKLJL @@ -78,7 +78,15 @@ Inputs and data of size N_samples x N_parameters (both arrays will be transposed kernel is used. The default kernel is the sum of a Squared Exponential kernel and white noise. """ -function GPObj(inputs, data, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} +function GPObj( + inputs, + data, + package::GPJL; + GPkernel::Union{K, KPy, Nothing}=nothing, + normalized::Bool=true, + prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} + + FT = eltype(data) models = Any[] @@ -117,11 +125,11 @@ function GPObj(inputs, data, package::GPJL; GPkernel::Union{K, KPy, Nothing}=not GPkernel_i = deepcopy(GPkernel) # inputs: N_samples x N_parameters # data: N_samples x N_data - logstd_obs_noise = log(sqrt(0.5)) # log standard dev of obs noise + logstd_regularization_noise = log(sqrt(0.5)) # log standard dev of obs noise # Zero mean function kmean = MeanZero() m = GPE(inputs', dropdims(data[:, i]', dims=1), kmean, GPkernel_i, - logstd_obs_noise) + logstd_regularization_noise) optimize!(m, noise=false) push!(models, m) end @@ -140,7 +148,13 @@ Inputs and data of size N_samples x N_parameters (both arrays will be transposed a default kernel is used. The default kernel is the sum of a Squared Exponential kernel and white noise. """ -function GPObj(inputs, data, package::SKLJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true) where {K<:Kernel, KPy<:PyObject} +function GPObj( + inputs, + data, + package::SKLJL; + GPkernel::Union{K, KPy, Nothing}=nothing, + normalized::Bool=true) where {K<:Kernel, KPy<:PyObject} + FT = eltype(data) models = Any[] @@ -186,6 +200,7 @@ function GPObj(inputs, data, package::SKLJL; GPkernel::Union{K, KPy, Nothing}=no end + """ predict(gp::GPObj{FT,GPJL}, new_inputs::Array{FT}) where {FT} = predict(gp, new_inputs, gp.prediction_type) @@ -234,4 +249,215 @@ function predict(gp::GPObj{FT, SKLJL}, return vcat(first.(μσ)...), vcat(last.(μσ)...).^2 end -end # module + +""" + NormalizedSVDGPObj{FT<:AbstractFloat} + +Structure holding training input and the fitted Gaussian Process Regression +models. + +# Fields +$(DocStringExtensions.FIELDS) + +""" +struct NormalizedSVDGPObj{FT<:AbstractFloat, GPM} + "training inputs/parameters; N_samples x N_parameters" + inputs::Array{FT} + "training data/targets; N_samples x N_data" + data::Array{FT} + "mean of input; 1 x N_parameters" + input_mean::Array{FT} + "square root of the inverse of the input covariance matrix; N_parameters x N_parameters" + sqrt_inv_input_cov::Union{Nothing, Array{FT, 2}} + "the Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs" + models::Vector + "whether to fit GP models on normalized inputs ((inputs - input_mean) * sqrt_inv_input_cov)" + normalized::Bool + "the singular value decomposition matrix" + decomposition + "the normalization in the svd transformed space" + sqrt_singular_values + "prediction type (`y` to predict the data, `f` to predict the latent function)" + prediction_type::Union{Nothing, PredictionType} +end + + +""" + GPObj(inputs, data, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, svdflag=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} + +Inputs and data of size N_samples x N_parameters (both arrays will be transposed in the construction of the GPObj) + + - `GPkernel` - GaussianProcesses kernel object. If not supplied, a default + kernel is used. The default kernel is the sum of a Squared + Exponential kernel and white noise. +""" + +function NormalizedSVDGPObj( + inputs, + data, + truth_cov::Array{AbstractFloat,2}, + package::GPJL; + GPkernel::Union{K, KPy, Nothing}=nothing, + normalized::Bool=true, + prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} + + # create an SVD decomposition of the covariance: + # NB: svdfact() may be more efficient / provide positive definiteness information + # stored as: svd.U * svd.S *svd.Vt + decomposition = svd(truth_cov) + + FT = eltype(data) + models = Any[] + + # Make sure inputs and data are arrays of type FT + inputs = convert(Array{FT}, inputs) + data = convert(Array{FT}, data) + + input_mean = reshape(mean(inputs, dims=1), 1, :) + sqrt_inv_input_cov = nothing + if normalized + if ndims(inputs) == 1 + error("Cov not defined for 1d input; can't normalize") + end + sqrt_inv_input_cov = convert(Array{FT}, sqrt(inv(cov(inputs, dims=1)))) + inputs = (inputs .- input_mean) * sqrt_inv_input_cov + end + + + sqrt_singular_values_inv=Diagonal(1.0 ./ sqrt.(SVD.S)) + transformed_data=sqrt_singular_values_inv * decomposition.Vt * data + + # Use a default kernel unless a kernel was supplied to GPObj + if GPkernel==nothing + # Construct kernel: + # Note that the kernels take the signal standard deviations on a + # log scale as input. + rbf_len = log.(ones(size(inputs, 2))) + rbf_logstd = log(1.0) + rbf = SEArd(rbf_len, rbf_logstd) + # regularize with white noise + white_logstd = log(1.0) + white = Noise(white_logstd) + # construct kernel + GPkernel = rbf + white + end + + for i in 1:size(transformed_data, 2) + # Make a copy of GPkernel (because the kernel gets altered in + # every iteration) + GPkernel_i = deepcopy(GPkernel) + # inputs: N_samples x N_parameters + # transformed_data: N_samples x N_data + + magic_number = 1e-3 # magic_number << 1 + logstd_regularization_noise = log(sqrt(magic_number)) # log standard dev of obs noise + + # Zero mean function + kmean = MeanZero() + m = GPE(inputs', + dropdims(transformed_data[:, i]', dims=1), + kmean, + GPkernel_i, + logstd_regularization_noise) + + # we choose above to explicitly learn the WhiteKernel as opposed to using this + # in built functionality. + optimize!(m, noise=false) + push!(models, m) + end + return NormalizedSVDGPObj{FT, typeof(package)}(inputs, + transformed_data, + input_mean, + sqrt_inv_input_cov, + models, + normalized, + decomposition, + inv(sqrt_singular_values_inv), + prediction_type) + +end + +""" + predict(gp::NormalizedSVDGPObj{FT,GPJL}, new_inputs::Array{FT}) where {FT} = predict(gp, new_inputs, gp.prediction_type) + +Evaluate the GP model(s) at new inputs. + - `gp` - a NormalizedSVDGPObj + - `new_inputs` - inputs for which GP model(s) is/are evaluated; N_new_inputs x N_parameters + +Note: If gp.normalized == true, the new inputs are normalized prior to the prediction +""" +predict( + gp::NormalizedSVDGPObj{FT, GPJL}, + new_inputs::Array{FT}; + transform_to_real=false) where {FT} = predict(gp, new_inputs, transform_to_real, gp.prediction_type) + +function predict( + gp::NormalizedSVDGPObj{FT, GPJL}, + new_inputs::Array{FT}, + transform_to_real, + ::FType) where {FT} + + if gp.normalized + new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov + end + M = length(gp.models) + μσ2 = [predict_f(gp.models[i], new_inputs) for i in 1:M] + # Return mean(s) and standard deviation(s) + μ = vcat(first.(μσ2)...) + σ2 = vcat(last.(μσ2)...) + if transform_to_real + #we created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP + μ = gp.decomposition.V * gp.sqrt_singular_values * μ + σ2 = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2) * gp.sqrt_singular_values * gp.decomposition.Vt + end + return μ, σ2 +end + +function predict( + gp::NormalizedSVDGPObj{FT, GPJL}, + new_inputs::Array{FT}, + transform_to_real, + ::YType) where {FT} + + if gp.normalized + new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov + end + M = length(gp.models) + # predicts columns of inputs so must be transposed + new_inputs = convert(Array{FT}, new_inputs') + μσ2 = [predict_y(gp.models[i], new_inputs) for i in 1:M] + # Return mean(s) and standard deviation(s) + μ = vcat(first.(μσ2)...) + σ2 = vcat(last.(μσ2)...) + if transform_to_real + #we created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP + μ = gp.decomposition.V * gp.sqrt_singular_values * μ + σ2 = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2) * gp.sqrt_singular_values * gp.decomposition.Vt + end + return μ, σ2 + +end + +function predict( + gp::NormalizedSVDGPObj{FT, SKLJL}, + new_inputs::Array{FT} + transform_to_real) where {FT} + + if gp.normalized + new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov + end + M = length(gp.models) + μσ = [gp.models[i].predict(new_inputs, return_std=true) for i in 1:M] + # Return mean(s) and standard deviation(s) + μ = vcat(first.(μσ2)...) + σ2 = vcat(last.(μσ2)...).^2 + if transform_to_real + #we created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP + μ = gp.decomposition.V * gp.sqrt_singular_values * μ + σ2 = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2) * gp.sqrt_singular_values * gp.decomposition.Vt + end + return μ, σ2 + +end + +end From a52a192c59eeff2a44fcae0b6deafcb6ba9cc7f5 Mon Sep 17 00:00:00 2001 From: odunbar Date: Mon, 29 Jun 2020 18:17:32 +0100 Subject: [PATCH 02/18] mcmc.jl now has svdflag, for transforming data sample into correct space --- src/MCMC.jl | 53 +++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 39 insertions(+), 14 deletions(-) diff --git a/src/MCMC.jl b/src/MCMC.jl index 79ef75863..e426aa10b 100644 --- a/src/MCMC.jl +++ b/src/MCMC.jl @@ -68,19 +68,34 @@ end where max_iter is the number of MCMC steps to perform (e.g., 100_000) """ -function MCMCObj(obs_sample::Vector{FT}, - obs_cov::Array{FT, 2}, - priors::Array{Prior, 1}, - step::FT, - param_init::Vector{FT}, - max_iter::IT, - algtype::String, - burnin::IT) where {FT<:AbstractFloat, IT<:Int} - +function MCMCObj( + obs_sample::Vector{FT}, + obs_cov::Array{FT, 2}, + priors::Array{Prior, 1}, + step::FT, + param_init::Vector{FT}, + max_iter::IT, + algtype::String, + burnin::IT; + svdflag=true) where {FT<:AbstractFloat, IT<:Int} + + + param_init_copy = deepcopy(param_init) + + #We need to transform obs_sample into the correct space ( + if svdflag + println("Applying SVD to decorrelating outputs, if not required set svdflag=false") + decomposition=svd(truth_cov)#svd.U * svd.S * svd.Vt (can also get V) + sqrt_singular_values_inv=Diagonal(1.0 ./ sqrt.(decomposition.S)) #diagonal matrix of 1/eigenvalues + obs_sample=sqrt_singular_values_inv*decomposition.Vt*obs_sample + else + println("Assuming independent outputs." + end + # first row is param_init - posterior = zeros(max_iter + 1, length(param_init)) - posterior[1, :] = param_init - param = param_init + posterior = zeros(max_iter + 1, length(param_init_copy)) + posterior[1, :] = param_init_copy + param = param_init_copy log_posterior = [nothing] iter = [1] obs_covinv = inv(obs_cov) @@ -88,8 +103,18 @@ function MCMCObj(obs_sample::Vector{FT}, if algtype != "rwm" error("only random walk metropolis 'rwm' is implemented so far") end - MCMCObj{FT,IT}(obs_sample, obs_cov, obs_covinv, priors, [step], burnin, - param, posterior, log_posterior, iter, accept, algtype) + MCMCObj{FT,IT}(obs_sample, + obs_cov, + obs_covinv, + priors, + [step], + burnin, + param, + posterior, + log_posterior, + iter, + accept, + algtype) end From 0fac07acb775eb8a21ea05211bb20cd1e83178aa Mon Sep 17 00:00:00 2001 From: odunbar Date: Mon, 29 Jun 2020 18:48:02 +0100 Subject: [PATCH 03/18] added noise-learn flag --- src/GPEmulator.jl | 72 +++++++++++++++++++++++++++++++++++++---------- 1 file changed, 57 insertions(+), 15 deletions(-) diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 12fa5de03..71539bfc2 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -84,8 +84,8 @@ function GPObj( package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, - prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} - + noise_learn=true, + prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} FT = eltype(data) models = Any[] @@ -106,35 +106,58 @@ function GPObj( # Use a default kernel unless a kernel was supplied to GPObj if GPkernel==nothing + println("Using default squared exponential kernel, learning lengthscale and variance parameters") # Construct kernel: # Note that the kernels take the signal standard deviations on a # log scale as input. rbf_len = log.(ones(size(inputs, 2))) rbf_logstd = log(1.0) rbf = SEArd(rbf_len, rbf_logstd) + kern=rbf + println("Using default squared exponential kernel:", kern) + + else + kern=deepcopy(GPkernel) + println("Using user-defined kernel",kern) + end + + + if noise_learn # regularize with white noise white_logstd = log(1.0) white = Noise(white_logstd) # construct kernel - GPkernel = rbf + white + kern = kern + white + println("Learning additive white noise") + + logstd_regularization_noise = log(sqrt(0.5)) # log standard dev of obs noise + else + logstd_regularization_noise = log(sqrt(1.0)) # log standard dev of obs noise end for i in 1:size(data, 2) # Make a copy of GPkernel (because the kernel gets altered in # every iteration) - GPkernel_i = deepcopy(GPkernel) + GPkernel_i = deepcopy(kern) # inputs: N_samples x N_parameters # data: N_samples x N_data - logstd_regularization_noise = log(sqrt(0.5)) # log standard dev of obs noise # Zero mean function kmean = MeanZero() - m = GPE(inputs', dropdims(data[:, i]', dims=1), kmean, GPkernel_i, + m = GPE(inputs', + dropdims(data[:, i]', dims=1), + kmean, + GPkernel_i, logstd_regularization_noise) + optimize!(m, noise=false) push!(models, m) end - return GPObj{FT, typeof(package)}(inputs, data, input_mean, - sqrt_inv_input_cov, models, normalized, + return GPObj{FT, typeof(package)}(inputs, + data, + input_mean, + sqrt_inv_input_cov, + models, + normalized, prediction_type) end @@ -299,6 +322,7 @@ function NormalizedSVDGPObj( package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, + noise_learn=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} # create an SVD decomposition of the covariance: @@ -329,29 +353,47 @@ function NormalizedSVDGPObj( # Use a default kernel unless a kernel was supplied to GPObj if GPkernel==nothing + println("Using default squared exponential kernel, learning lengthscale and variance parameters") # Construct kernel: # Note that the kernels take the signal standard deviations on a # log scale as input. rbf_len = log.(ones(size(inputs, 2))) rbf_logstd = log(1.0) rbf = SEArd(rbf_len, rbf_logstd) + + kern=rbf + println("Using default squared exponential kernel:", kern) + else + kern=deepcopy(GPkernel) + println("Using user-defined kernel",kern) + end + + + if noise_learn # regularize with white noise white_logstd = log(1.0) white = Noise(white_logstd) - # construct kernel - GPkernel = rbf + white - end + # add to kernel construct kernel + kern = kern + white + println("Learning additive white noise") + + #make the regularization small. We actually learn + # total_noise = white_logstd + logstd_regularization_noise) + magic_number = 1e-3 # magic_number << 1 + logstd_regularization_noise = log(sqrt(1e-3)) + else + # when not learning noise, our SVD transformation implies the observational noise is the identity. + logstd_regularization_noise = log(sqrt(1.0)) + end + for i in 1:size(transformed_data, 2) # Make a copy of GPkernel (because the kernel gets altered in # every iteration) - GPkernel_i = deepcopy(GPkernel) + GPkernel_i = deepcopy(kern) # inputs: N_samples x N_parameters # transformed_data: N_samples x N_data - magic_number = 1e-3 # magic_number << 1 - logstd_regularization_noise = log(sqrt(magic_number)) # log standard dev of obs noise - # Zero mean function kmean = MeanZero() m = GPE(inputs', From c4e6e2e56dd125e995e0cbfda5da2c708a910d6d Mon Sep 17 00:00:00 2001 From: odunbar Date: Mon, 29 Jun 2020 19:03:09 +0100 Subject: [PATCH 04/18] replaced old GPObj, updated the SKlearn constructor --- src/GPEmulator.jl | 306 ++++++++++++++-------------------------------- 1 file changed, 92 insertions(+), 214 deletions(-) diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 71539bfc2..627b9181e 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -64,13 +64,17 @@ struct GPObj{FT<:AbstractFloat, GPM} models::Vector "whether to fit GP models on normalized inputs ((inputs - input_mean) * sqrt_inv_input_cov)" normalized::Bool + "the singular value decomposition matrix" + decomposition + "the normalization in the svd transformed space" + sqrt_singular_values "prediction type (`y` to predict the data, `f` to predict the latent function)" prediction_type::Union{Nothing, PredictionType} end """ - GPObj(inputs, data, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} + GPObj(inputs, data, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, svdflag=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} Inputs and data of size N_samples x N_parameters (both arrays will be transposed in the construction of the GPObj) @@ -78,22 +82,31 @@ Inputs and data of size N_samples x N_parameters (both arrays will be transposed kernel is used. The default kernel is the sum of a Squared Exponential kernel and white noise. """ + function GPObj( inputs, data, + truth_cov::Array{AbstractFloat,2}, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, noise_learn=true, - prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} + prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} + FT = eltype(data) models = Any[] + # create an SVD decomposition of the covariance: + # NB: svdfact() may be more efficient / provide positive definiteness information + # stored as: svd.U * svd.S *svd.Vt + decomposition = svd(truth_cov) + # Make sure inputs and data are arrays of type FT inputs = convert(Array{FT}, inputs) data = convert(Array{FT}, data) + # Transform the inputs input_mean = reshape(mean(inputs, dims=1), 1, :) sqrt_inv_input_cov = nothing if normalized @@ -103,7 +116,11 @@ function GPObj( sqrt_inv_input_cov = convert(Array{FT}, sqrt(inv(cov(inputs, dims=1)))) inputs = (inputs .- input_mean) * sqrt_inv_input_cov end - + + # Transform the outputs + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(SVD.S)) + transformed_data = sqrt_singular_values_inv * decomposition.Vt * data + # Use a default kernel unless a kernel was supplied to GPObj if GPkernel==nothing println("Using default squared exponential kernel, learning lengthscale and variance parameters") @@ -113,9 +130,9 @@ function GPObj( rbf_len = log.(ones(size(inputs, 2))) rbf_logstd = log(1.0) rbf = SEArd(rbf_len, rbf_logstd) + kern=rbf println("Using default squared exponential kernel:", kern) - else kern=deepcopy(GPkernel) println("Using user-defined kernel",kern) @@ -126,42 +143,54 @@ function GPObj( # regularize with white noise white_logstd = log(1.0) white = Noise(white_logstd) - # construct kernel + + # add to kernel construct kernel kern = kern + white println("Learning additive white noise") - logstd_regularization_noise = log(sqrt(0.5)) # log standard dev of obs noise + #make the regularization small. We actually learn + # total_noise = white_logstd + logstd_regularization_noise) + magic_number = 1e-3 # magic_number << 1 + logstd_regularization_noise = log(sqrt(1e-3)) else - logstd_regularization_noise = log(sqrt(1.0)) # log standard dev of obs noise + # when not learning noise, our SVD transformation implies the observational noise is the identity. + logstd_regularization_noise = log(sqrt(1.0)) end - - for i in 1:size(data, 2) + + for i in 1:size(transformed_data, 2) # Make a copy of GPkernel (because the kernel gets altered in # every iteration) GPkernel_i = deepcopy(kern) - # inputs: N_samples x N_parameters - # data: N_samples x N_data + # inputs: N_samples x N_parameters + # transformed_data: N_samples x N_data + # Zero mean function kmean = MeanZero() m = GPE(inputs', - dropdims(data[:, i]', dims=1), + dropdims(transformed_data[:, i]', dims=1), kmean, GPkernel_i, logstd_regularization_noise) - + + # we choose above to explicitly learn the WhiteKernel as opposed to using this + # in built functionality. optimize!(m, noise=false) push!(models, m) end return GPObj{FT, typeof(package)}(inputs, - data, + transformed_data, input_mean, sqrt_inv_input_cov, models, normalized, + decomposition, + inv(sqrt_singular_values_inv), prediction_type) + end + """ GPObj(inputs, data, package::SKLJL, GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true) where {K<:Kernel, KPy<:PyObject} @@ -174,6 +203,7 @@ Inputs and data of size N_samples x N_parameters (both arrays will be transposed function GPObj( inputs, data, + truth_cov::Array{AbstractFloat,2}, package::SKLJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true) where {K<:Kernel, KPy<:PyObject} @@ -181,162 +211,16 @@ function GPObj( FT = eltype(data) models = Any[] - # Make sure inputs and data are arrays of type FT - inputs = convert(Array{FT}, inputs) - data = convert(Array{FT}, data) - - input_mean = reshape(mean(inputs, dims=1), 1, :) - sqrt_inv_input_cov = nothing - if normalized - if ndims(inputs) == 1 - error("Cov not defined for 1d input; can't normalize") - end - sqrt_inv_input_cov = convert(Array{FT}, sqrt(inv(cov(inputs, dims=1)))) - inputs = (inputs .- input_mean) * sqrt_inv_input_cov - end - - if GPkernel==nothing - const_value = 1.0 - var_kern = ConstantKernel(constant_value=const_value, - constant_value_bounds=(1e-05, 10000.0)) - rbf_len = ones(size(inputs, 2)) - rbf = RBF(length_scale=rbf_len, length_scale_bounds=(1e-05, 10000.0)) - white_noise_level = 1.0 - white = WhiteKernel(noise_level=white_noise_level, - noise_level_bounds=(1e-05, 10.0)) - GPkernel = var_kern * rbf + white - end - for i in 1:size(data, 2) - m = GaussianProcessRegressor(kernel=GPkernel, - n_restarts_optimizer=10, - alpha=0.0, normalize_y=true) - ScikitLearn.fit!(m, inputs, data[:, i]) - if i==1 - println(m.kernel.hyperparameters) - print("Completed training of: ") - end - print(i,", ") - push!(models, m) - end - return GPObj{FT, typeof(package)}(inputs, data, input_mean, sqrt_inv_input_cov, - models, normalized, YType()) -end - - - -""" - predict(gp::GPObj{FT,GPJL}, new_inputs::Array{FT}) where {FT} = predict(gp, new_inputs, gp.prediction_type) - -Evaluate the GP model(s) at new inputs. - - `gp` - a GPObj - - `new_inputs` - inputs for which GP model(s) is/are evaluated; N_new_inputs x N_parameters - -Note: If gp.normalized == true, the new inputs are normalized prior to the prediction -""" -predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}) where {FT} = predict(gp, new_inputs, gp.prediction_type) - -function predict(gp::GPObj{FT, GPJL}, - new_inputs::Array{FT}, - ::FType) where {FT} - if gp.normalized - new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov - end - M = length(gp.models) - μσ2 = [predict_f(gp.models[i], new_inputs) for i in 1:M] - # Return mean(s) and standard deviation(s) - return vcat(first.(μσ2)...), vcat(last.(μσ2)...) -end - -function predict(gp::GPObj{FT, GPJL}, - new_inputs::Array{FT}, - ::YType) where {FT} - if gp.normalized - new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov - end - M = length(gp.models) - # predicts columns of inputs so must be transposed - new_inputs = convert(Array{FT}, new_inputs') - μσ2 = [predict_y(gp.models[i], new_inputs) for i in 1:M] - # Return mean(s) and standard deviation(s) - return vcat(first.(μσ2)...), vcat(last.(μσ2)...) -end - -function predict(gp::GPObj{FT, SKLJL}, - new_inputs::Array{FT}) where {FT} - if gp.normalized - new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov - end - M = length(gp.models) - μσ = [gp.models[i].predict(new_inputs, return_std=true) for i in 1:M] - # Return mean(s) and standard deviation(s) - return vcat(first.(μσ)...), vcat(last.(μσ)...).^2 -end - - -""" - NormalizedSVDGPObj{FT<:AbstractFloat} - -Structure holding training input and the fitted Gaussian Process Regression -models. - -# Fields -$(DocStringExtensions.FIELDS) - -""" -struct NormalizedSVDGPObj{FT<:AbstractFloat, GPM} - "training inputs/parameters; N_samples x N_parameters" - inputs::Array{FT} - "training data/targets; N_samples x N_data" - data::Array{FT} - "mean of input; 1 x N_parameters" - input_mean::Array{FT} - "square root of the inverse of the input covariance matrix; N_parameters x N_parameters" - sqrt_inv_input_cov::Union{Nothing, Array{FT, 2}} - "the Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs" - models::Vector - "whether to fit GP models on normalized inputs ((inputs - input_mean) * sqrt_inv_input_cov)" - normalized::Bool - "the singular value decomposition matrix" - decomposition - "the normalization in the svd transformed space" - sqrt_singular_values - "prediction type (`y` to predict the data, `f` to predict the latent function)" - prediction_type::Union{Nothing, PredictionType} -end - - -""" - GPObj(inputs, data, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, svdflag=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} - -Inputs and data of size N_samples x N_parameters (both arrays will be transposed in the construction of the GPObj) - - - `GPkernel` - GaussianProcesses kernel object. If not supplied, a default - kernel is used. The default kernel is the sum of a Squared - Exponential kernel and white noise. -""" - -function NormalizedSVDGPObj( - inputs, - data, - truth_cov::Array{AbstractFloat,2}, - package::GPJL; - GPkernel::Union{K, KPy, Nothing}=nothing, - normalized::Bool=true, - noise_learn=true, - prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} - # create an SVD decomposition of the covariance: # NB: svdfact() may be more efficient / provide positive definiteness information # stored as: svd.U * svd.S *svd.Vt decomposition = svd(truth_cov) - - FT = eltype(data) - models = Any[] - + # Make sure inputs and data are arrays of type FT inputs = convert(Array{FT}, inputs) data = convert(Array{FT}, data) + # Transform the inputs input_mean = reshape(mean(inputs, dims=1), 1, :) sqrt_inv_input_cov = nothing if normalized @@ -347,94 +231,88 @@ function NormalizedSVDGPObj( inputs = (inputs .- input_mean) * sqrt_inv_input_cov end + # Transform the outputs + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(SVD.S)) + transformed_data = sqrt_singular_values_inv * decomposition.Vt * data - sqrt_singular_values_inv=Diagonal(1.0 ./ sqrt.(SVD.S)) - transformed_data=sqrt_singular_values_inv * decomposition.Vt * data - - # Use a default kernel unless a kernel was supplied to GPObj if GPkernel==nothing println("Using default squared exponential kernel, learning lengthscale and variance parameters") - # Construct kernel: - # Note that the kernels take the signal standard deviations on a - # log scale as input. - rbf_len = log.(ones(size(inputs, 2))) - rbf_logstd = log(1.0) - rbf = SEArd(rbf_len, rbf_logstd) + #Create default squared exponential kernel + const_value = 1.0 + var_kern = ConstantKernel(constant_value=const_value, + constant_value_bounds=(1e-05, 10000.0)) + rbf_len = ones(size(inputs, 2)) + rbf = RBF(length_scale=rbf_len, length_scale_bounds=(1e-05, 10000.0)) - kern=rbf + + kern = var_kern * rbf println("Using default squared exponential kernel:", kern) else - kern=deepcopy(GPkernel) + kern = deepcopy(GPkernel) println("Using user-defined kernel",kern) end - if noise_learn # regularize with white noise - white_logstd = log(1.0) - white = Noise(white_logstd) + white_noise_level = 1.0 + white = WhiteKernel(noise_level=white_noise_level, + noise_level_bounds=(1e-05, 10.0)) - # add to kernel construct kernel kern = kern + white println("Learning additive white noise") #make the regularization small. We actually learn - # total_noise = white_logstd + logstd_regularization_noise) + # total_noise = white_noise_level + regularization_noise magic_number = 1e-3 # magic_number << 1 - logstd_regularization_noise = log(sqrt(1e-3)) + regularization_noise = 1e-3 else # when not learning noise, our SVD transformation implies the observational noise is the identity. - logstd_regularization_noise = log(sqrt(1.0)) + regularization_noise = 1.0 end + - for i in 1:size(transformed_data, 2) - # Make a copy of GPkernel (because the kernel gets altered in - # every iteration) + for i in 1:size(data, 2) GPkernel_i = deepcopy(kern) - # inputs: N_samples x N_parameters - # transformed_data: N_samples x N_data - - # Zero mean function - kmean = MeanZero() - m = GPE(inputs', - dropdims(transformed_data[:, i]', dims=1), - kmean, - GPkernel_i, - logstd_regularization_noise) - - # we choose above to explicitly learn the WhiteKernel as opposed to using this - # in built functionality. - optimize!(m, noise=false) + + m = GaussianProcessRegressor(kernel=GPkernel_i, + n_restarts_optimizer=10, + alpha=regularization_noise) + + ScikitLearn.fit!(m, inputs, data[:, i]) + if i==1 + println(m.kernel.hyperparameters) + print("Completed training of: ") + end + print(i,", ") push!(models, m) end - return NormalizedSVDGPObj{FT, typeof(package)}(inputs, - transformed_data, - input_mean, - sqrt_inv_input_cov, - models, - normalized, - decomposition, - inv(sqrt_singular_values_inv), - prediction_type) - + return GPObj{FT, typeof(package)}(inputs, + data, + input_mean, + sqrt_inv_input_cov, + models, + normalized, + YType()) end + + """ - predict(gp::NormalizedSVDGPObj{FT,GPJL}, new_inputs::Array{FT}) where {FT} = predict(gp, new_inputs, gp.prediction_type) + predict(gp::GPObj{FT,GPJL}, new_inputs::Array{FT}) where {FT} = predict(gp, new_inputs, gp.prediction_type) Evaluate the GP model(s) at new inputs. - - `gp` - a NormalizedSVDGPObj + - `gp` - a GPObj - `new_inputs` - inputs for which GP model(s) is/are evaluated; N_new_inputs x N_parameters Note: If gp.normalized == true, the new inputs are normalized prior to the prediction """ predict( - gp::NormalizedSVDGPObj{FT, GPJL}, + gp::GPObj{FT, GPJL}, new_inputs::Array{FT}; transform_to_real=false) where {FT} = predict(gp, new_inputs, transform_to_real, gp.prediction_type) function predict( - gp::NormalizedSVDGPObj{FT, GPJL}, + gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, ::FType) where {FT} @@ -456,7 +334,7 @@ function predict( end function predict( - gp::NormalizedSVDGPObj{FT, GPJL}, + gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, ::YType) where {FT} @@ -481,7 +359,7 @@ function predict( end function predict( - gp::NormalizedSVDGPObj{FT, SKLJL}, + gp::GPObj{FT, SKLJL}, new_inputs::Array{FT} transform_to_real) where {FT} From 68d6c78ca7a4d2a244d5cc44de5e2ce6415ed181 Mon Sep 17 00:00:00 2001 From: odunbar Date: Mon, 29 Jun 2020 20:02:17 +0100 Subject: [PATCH 05/18] typos --- src/GPEmulator.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 627b9181e..542cabb3e 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -65,9 +65,9 @@ struct GPObj{FT<:AbstractFloat, GPM} "whether to fit GP models on normalized inputs ((inputs - input_mean) * sqrt_inv_input_cov)" normalized::Bool "the singular value decomposition matrix" - decomposition - "the normalization in the svd transformed space" - sqrt_singular_values + decomposition::SVD + "the normalization matrix in the svd transformed space" + sqrt_singular_values::Array{FT,2} "prediction type (`y` to predict the data, `f` to predict the latent function)" prediction_type::Union{Nothing, PredictionType} end @@ -118,7 +118,7 @@ function GPObj( end # Transform the outputs - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(SVD.S)) + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) transformed_data = sqrt_singular_values_inv * decomposition.Vt * data # Use a default kernel unless a kernel was supplied to GPObj @@ -232,7 +232,7 @@ function GPObj( end # Transform the outputs - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(SVD.S)) + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) transformed_data = sqrt_singular_values_inv * decomposition.Vt * data if GPkernel==nothing From cb397ca600decd3bd9525890abc3216dd4518573 Mon Sep 17 00:00:00 2001 From: Melanie Bieli Date: Wed, 1 Jul 2020 19:45:56 -0400 Subject: [PATCH 06/18] Add example of a GP plot --- examples/GPEmulator/plot_GP.jl | 94 ++++++++++++++++++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 examples/GPEmulator/plot_GP.jl diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl new file mode 100644 index 000000000..957bbdbb3 --- /dev/null +++ b/examples/GPEmulator/plot_GP.jl @@ -0,0 +1,94 @@ +# Import modules +using Random +using GaussianProcesses +using Distributions +using Statistics +using Plots; pyplot(size = (1000, 400)) + +using CalibrateEmulateSample.GPEmulator + +############################################################################### +# # +# This examples shows how to fit a Gaussian Process regression model to # +# 2D input data, and how to plot the mean and variance # +# # +############################################################################### + +function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T + m, n = length(vy), length(vx) + gx = reshape(repeat(vx, inner = m, outer = 1), m, n) + gy = reshape(repeat(vy, inner = 1, outer = n), m, n) + + return gx, gy +end + +gppackage = GPJL() +pred_type = YType() + +# Seed for pseudo-random number generator +rng_seed = 41 +Random.seed!(rng_seed) + + +# Generate training data (x-y pairs, where x ∈ ℝ ᵖ, y ∈ ℝ ᵈ) +# x = [x1, x2]: inputs/predictors/features/parameters +# y = [y1, y2]: outputs/predictands/targets/observables +# The observables y are related to the parameters x by: +# y = G(x1, x2) + η, +# where G(x1, x2) := [sin(x1) + cos(x2), sin(x1) - cos(x2)], and η ~ N(0, Σ) +n = 100 # number of training points +p = 2 # input dim +d = 2 # output dim +X = 2.0 * π * rand(n, p) +μ = zeros(d) +Σ = 0.1 * [[0.5, 0.2] [0.2, 0.8]] # d x d + +g1 = sin.(X[:, 1]) .+ cos.(X[:, 2]) +g2 = sin.(X[:, 1]) .- cos.(X[:, 2]) +Y = [g1 g2] .+ rand(MvNormal(μ, Σ), n)' +truth_cov = cov(Y, dims=1) + +# Fit 2D Gaussian Process regression model +# (To be precise, we fit two models, one that predicts y1 from x1 and x2, +# and one that predicts y2 from x1 and x2) +gpobj = GPObj(X, Y, truth_cov, gppackage, GPkernel=nothing, + normalized=true, noise_learn=true, prediction_type=pred_type) + +# Plot mean and variance of the predicted observables y1 and y2 +# For this, we generate test points on a x1-x2 grid +n_pts = 50 +x1 = range(0.0, stop=2*π, length=n_pts) +x2 = range(0.0, stop=2*π, length=n_pts) +X1, X2 = meshgrid(x1, x2) + +# Input for predict needs to be N_samples x N_parameters +inputs = hcat(X1[:], X2[:]) + +# Get predictions on the grid +gp_mean, gp_var = GPEmulator.predict(gpobj, inputs, transform_to_real=true) + +for y_i in 1:d + mean_grid = reshape(gp_mean[y_i, :], n_pts, n_pts) + p1 = plot(x1, x2, mean_grid, st=:surface, camera=(-30, 30), c=:cividis, + xlabel="x1", ylabel="x2", zlabel="mean of y"*string(y_i)) + + var_grid = reshape(gp_var[y_i, :], n_pts, n_pts) + p2 = plot(x1, x2, var_grid, st=:surface, camera=(-30, 30), c=:cividis, + xlabel="x1", ylabel="x2", zlabel="var of y"*string(y_i)) + + plot(p1, p2, layout=(1, 2), legend=false) + savefig("/home/melanie/Desktop/GP_test_y"*string(y_i)*".png") +end + +# Make plots of the true components of G(x1, x2) +g1_true = sin.(inputs[:, 1]) .+ cos.(inputs[:, 2]) +g1_grid = reshape(g1_true, n_pts, n_pts) +p3 = plot(x1, x2, g1_grid, st=:surface, camera=(-30, 30), c=:cividis, + xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") +savefig("/home/melanie/Desktop/GP_test_true_g1.png") + +g2_true = sin.(inputs[:, 1]) .- cos.(inputs[:, 2]) +g2_grid = reshape(g2_true, n_pts, n_pts) +p4 = plot(x1, x2, g2_grid, st=:surface, camera=(-30, 30), c=:cividis, + xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") +savefig("/home/melanie/Desktop/GP_test_true_g2.png") From 98f707a0f49402c555c5041b7946196866cdc155 Mon Sep 17 00:00:00 2001 From: bielim Date: Thu, 2 Jul 2020 07:01:49 -0700 Subject: [PATCH 07/18] Update examples/GPEmulator/plot_GP.jl Change path where figure gets saved Co-authored-by: Charles Kawczynski --- examples/GPEmulator/plot_GP.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl index 957bbdbb3..2493d5734 100644 --- a/examples/GPEmulator/plot_GP.jl +++ b/examples/GPEmulator/plot_GP.jl @@ -85,7 +85,7 @@ g1_true = sin.(inputs[:, 1]) .+ cos.(inputs[:, 2]) g1_grid = reshape(g1_true, n_pts, n_pts) p3 = plot(x1, x2, g1_grid, st=:surface, camera=(-30, 30), c=:cividis, xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") -savefig("/home/melanie/Desktop/GP_test_true_g1.png") +savefig("GP_test_true_g1.png") g2_true = sin.(inputs[:, 1]) .- cos.(inputs[:, 2]) g2_grid = reshape(g2_true, n_pts, n_pts) From 069ab5436371b57f99dbb9ba217562f7c58466ef Mon Sep 17 00:00:00 2001 From: bielim Date: Thu, 2 Jul 2020 07:02:48 -0700 Subject: [PATCH 08/18] Update examples/GPEmulator/plot_GP.jl Change path where figure gets saved Co-authored-by: Charles Kawczynski --- examples/GPEmulator/plot_GP.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl index 2493d5734..7d2a0cab4 100644 --- a/examples/GPEmulator/plot_GP.jl +++ b/examples/GPEmulator/plot_GP.jl @@ -91,4 +91,4 @@ g2_true = sin.(inputs[:, 1]) .- cos.(inputs[:, 2]) g2_grid = reshape(g2_true, n_pts, n_pts) p4 = plot(x1, x2, g2_grid, st=:surface, camera=(-30, 30), c=:cividis, xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") -savefig("/home/melanie/Desktop/GP_test_true_g2.png") +savefig("GP_test_true_g2.png") From 1661de0d8e7f61e8eb54db4b1705e54014f05f5d Mon Sep 17 00:00:00 2001 From: Melanie Bieli Date: Thu, 2 Jul 2020 10:44:48 -0400 Subject: [PATCH 09/18] Try to make the `predict` function work for matrices of test points --- src/GPEmulator.jl | 126 +++++++++++++++++++++++++++++----------------- 1 file changed, 80 insertions(+), 46 deletions(-) diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 542cabb3e..3bf27d722 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -84,17 +84,17 @@ Inputs and data of size N_samples x N_parameters (both arrays will be transposed """ function GPObj( - inputs, - data, - truth_cov::Array{AbstractFloat,2}, + inputs::Array{FT, 2}, + data::Array{FT, 2}, + truth_cov::Array{FT, 2}, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, - noise_learn=true, - prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} + noise_learn::Bool=true, + prediction_type::PredictionType=YType()) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject} - FT = eltype(data) + #FT = eltype(data) models = Any[] # create an SVD decomposition of the covariance: @@ -119,7 +119,7 @@ function GPObj( # Transform the outputs sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - transformed_data = sqrt_singular_values_inv * decomposition.Vt * data + transformed_data = Array((sqrt_singular_values_inv * decomposition.Vt * data')') # Use a default kernel unless a kernel was supplied to GPObj if GPkernel==nothing @@ -138,7 +138,6 @@ function GPObj( println("Using user-defined kernel",kern) end - if noise_learn # regularize with white noise white_logstd = log(1.0) @@ -177,6 +176,7 @@ function GPObj( optimize!(m, noise=false) push!(models, m) end + return GPObj{FT, typeof(package)}(inputs, transformed_data, input_mean, @@ -186,7 +186,6 @@ function GPObj( decomposition, inv(sqrt_singular_values_inv), prediction_type) - end @@ -201,14 +200,16 @@ Inputs and data of size N_samples x N_parameters (both arrays will be transposed Squared Exponential kernel and white noise. """ function GPObj( - inputs, - data, - truth_cov::Array{AbstractFloat,2}, + inputs::Array{FT, 2}, + data::Array{FT, 2}, + truth_cov::Array{FT, 2}, package::SKLJL; GPkernel::Union{K, KPy, Nothing}=nothing, - normalized::Bool=true) where {K<:Kernel, KPy<:PyObject} + normalized::Bool=true, + noise_learn::Bool=true, + prediction_type::PredictionType=YType()) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject} - FT = eltype(data) + #FT = eltype(data) models = Any[] # create an SVD decomposition of the covariance: @@ -233,7 +234,8 @@ function GPObj( # Transform the outputs sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - transformed_data = sqrt_singular_values_inv * decomposition.Vt * data + #transformed_data = sqrt_singular_values_inv * decomposition.Vt * data + transformed_data = Array((sqrt_singular_values_inv * decomposition.Vt * data')') if GPkernel==nothing println("Using default squared exponential kernel, learning lengthscale and variance parameters") @@ -292,6 +294,8 @@ function GPObj( sqrt_inv_input_cov, models, normalized, + decomposition, + inv(sqrt_singular_values_inv), YType()) end @@ -306,16 +310,10 @@ Evaluate the GP model(s) at new inputs. Note: If gp.normalized == true, the new inputs are normalized prior to the prediction """ -predict( - gp::GPObj{FT, GPJL}, - new_inputs::Array{FT}; - transform_to_real=false) where {FT} = predict(gp, new_inputs, transform_to_real, gp.prediction_type) +predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}; transform_to_real=false) where {FT} = predict(gp, new_inputs, transform_to_real, gp.prediction_type) -function predict( - gp::GPObj{FT, GPJL}, - new_inputs::Array{FT}, - transform_to_real, - ::FType) where {FT} +function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, + ::FType) where {FT} if gp.normalized new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov @@ -323,22 +321,31 @@ function predict( M = length(gp.models) μσ2 = [predict_f(gp.models[i], new_inputs) for i in 1:M] # Return mean(s) and standard deviation(s) - μ = vcat(first.(μσ2)...) - σ2 = vcat(last.(μσ2)...) + μ = reshape(vcat(first.(μσ2)...), + size(new_inputs)[1], size(new_inputs)[2]) + σ2 = reshape(vcat(last.(μσ2)...), + size(new_inputs)[1], size(new_inputs)[2]) if transform_to_real - #we created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP + # Revert back to the original (correlated) coordinate system + # We created meanvGP = Dinv*Vt*meanv, so meanv = V*D*meanvGP μ = gp.decomposition.V * gp.sqrt_singular_values * μ - σ2 = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2) * gp.sqrt_singular_values * gp.decomposition.Vt + temp = zeros(size(new_inputs)) + # Back transformation of covariance + for j in 1:size(new_inputs, 2) + σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[:,j]) * gp.sqrt_singular_values * gp.decomposition.Vt + # Extract diagonal elements from σ2_j and add them to σ2 as column + # vectors + # TODO: Should we have an optional input argument `full_cov` - if + # set to true, the full covariance matrix is returned instead of + # only variances (default: false)? + temp[:, j] = diag(σ2_j) + end + σ2 = temp end return μ, σ2 end -function predict( - gp::GPObj{FT, GPJL}, - new_inputs::Array{FT}, - transform_to_real, - ::YType) where {FT} - +function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, ::YType) where {FT} if gp.normalized new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov end @@ -346,22 +353,34 @@ function predict( # predicts columns of inputs so must be transposed new_inputs = convert(Array{FT}, new_inputs') μσ2 = [predict_y(gp.models[i], new_inputs) for i in 1:M] + # Return mean(s) and standard deviation(s) - μ = vcat(first.(μσ2)...) - σ2 = vcat(last.(μσ2)...) + μ = reshape(vcat(first.(μσ2)...), + size(new_inputs)[1], size(new_inputs)[2]) + σ2 = reshape(vcat(last.(μσ2)...), + size(new_inputs)[1], size(new_inputs)[2]) if transform_to_real - #we created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP + # Revert back to the original (correlated) coordinate system + # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP μ = gp.decomposition.V * gp.sqrt_singular_values * μ - σ2 = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2) * gp.sqrt_singular_values * gp.decomposition.Vt + temp = zeros(size(σ2)) + # Back transformation of covariance + for j in 1:size(σ2, 2) + σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[:, j]) * gp.sqrt_singular_values * gp.decomposition.Vt + # Extract diagonal elements from σ2_j and add them to σ2 as column + # vectors + # TODO: Should we have an optional input argument `full_cov` - if + # set to true, the full covariance matrix is returned instead of + # only variances (default: false)? + temp[:, j] = diag(σ2_j) + end + σ2 = temp end return μ, σ2 end -function predict( - gp::GPObj{FT, SKLJL}, - new_inputs::Array{FT} - transform_to_real) where {FT} +function predict(gp::GPObj{FT, SKLJL}, new_inputs::Array{FT}, transform_to_real) where {FT} if gp.normalized new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov @@ -369,12 +388,27 @@ function predict( M = length(gp.models) μσ = [gp.models[i].predict(new_inputs, return_std=true) for i in 1:M] # Return mean(s) and standard deviation(s) - μ = vcat(first.(μσ2)...) - σ2 = vcat(last.(μσ2)...).^2 + μ = reshape(vcat(first.(μσ2)...), + size(new_inputs)[1], size(new_inputs)[2]) + σ2 = reshape(vcat(last.(μσ2)...), + size(new_inputs)[1], size(new_inputs)[2]) + if transform_to_real - #we created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP + # Revert back to the original (correlated) coordinate system + # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP μ = gp.decomposition.V * gp.sqrt_singular_values * μ - σ2 = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2) * gp.sqrt_singular_values * gp.decomposition.Vt + temp = zeros(size(σ2)) + # Back transformation of covariance + for j in 1:size(σ2, 2) + σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[:, j]) * gp.sqrt_singular_values * gp.decomposition.Vt + # Extract diagonal elements from σ2_j and add them to σ2 as column + # vectors + # TODO: Should we have an optional input argument `full_cov` - if + # set to true, the full covariance matrix is returned instead of + # only variances (default: false)? + temp[:, j] = diag(σ2_j) + end + σ2 = temp end return μ, σ2 From 4c2cee83d50f84f6f8508d26116100115122627c Mon Sep 17 00:00:00 2001 From: odunbar Date: Tue, 7 Jul 2020 17:57:10 +0100 Subject: [PATCH 10/18] example improvements, bugfixes --- examples/GPEmulator/plot_GP.jl | 29 +++++++++++----- src/GPEmulator.jl | 63 +++++++++++++++++++++------------- src/MCMC.jl | 2 +- 3 files changed, 60 insertions(+), 34 deletions(-) diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl index 7d2a0cab4..2302d0c7b 100644 --- a/examples/GPEmulator/plot_GP.jl +++ b/examples/GPEmulator/plot_GP.jl @@ -4,7 +4,7 @@ using GaussianProcesses using Distributions using Statistics using Plots; pyplot(size = (1000, 400)) - +using LinearAlgebra using CalibrateEmulateSample.GPEmulator ############################################################################### @@ -39,20 +39,31 @@ Random.seed!(rng_seed) n = 100 # number of training points p = 2 # input dim d = 2 # output dim -X = 2.0 * π * rand(n, p) +#X = 2.0 * π * rand(n, p) +x1 = range(0.0, stop=2*π, length=convert(Int64,sqrt(n))) +x2 = range(0.0, stop=2*π, length=convert(Int64,sqrt(n))) +X1, X2 = meshgrid(x1, x2) + +# Input for predict needs to be N_samples x N_parameters +X = hcat(X1[:], X2[:]) μ = zeros(d) -Σ = 0.1 * [[0.5, 0.2] [0.2, 0.8]] # d x d +Σ = 0.1 * [[0.8, 0.0] [0.0, 0.5]] # d x d + +g1x = sin.(X[:, 1]) .+ cos.(X[:, 2]) +g2x = sin.(X[:, 1]) .- cos.(X[:, 2]) +noise_samples = rand(MvNormal(μ, Σ), n)' +gx = [g1x g2x] +Y = gx .+ noise_samples + +truth_cov = Σ +sample_truth_cov = cov(noise_samples, dims=1) -g1 = sin.(X[:, 1]) .+ cos.(X[:, 2]) -g2 = sin.(X[:, 1]) .- cos.(X[:, 2]) -Y = [g1 g2] .+ rand(MvNormal(μ, Σ), n)' -truth_cov = cov(Y, dims=1) # Fit 2D Gaussian Process regression model # (To be precise, we fit two models, one that predicts y1 from x1 and x2, # and one that predicts y2 from x1 and x2) gpobj = GPObj(X, Y, truth_cov, gppackage, GPkernel=nothing, - normalized=true, noise_learn=true, prediction_type=pred_type) + normalized=false, noise_learn=false, prediction_type=pred_type) # Plot mean and variance of the predicted observables y1 and y2 # For this, we generate test points on a x1-x2 grid @@ -77,7 +88,7 @@ for y_i in 1:d xlabel="x1", ylabel="x2", zlabel="var of y"*string(y_i)) plot(p1, p2, layout=(1, 2), legend=false) - savefig("/home/melanie/Desktop/GP_test_y"*string(y_i)*".png") + savefig("GP_test_y"*string(y_i)*".png") end # Make plots of the true components of G(x1, x2) diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 3bf27d722..355ae0400 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -12,7 +12,7 @@ using PyCall @sk_import gaussian_process : GaussianProcessRegressor @sk_import gaussian_process.kernels : (RBF, WhiteKernel, ConstantKernel) -export GPObj, NormalizedSVDGPObj +export GPObj export predict export GPJL, SKLJL @@ -114,10 +114,13 @@ function GPObj( error("Cov not defined for 1d input; can't normalize") end sqrt_inv_input_cov = convert(Array{FT}, sqrt(inv(cov(inputs, dims=1)))) - inputs = (inputs .- input_mean) * sqrt_inv_input_cov + new_inputs = (inputs .- input_mean) * sqrt_inv_input_cov + else + new_inputs = inputs end - + # Transform the outputs + #data = convert(Matrix{Float64},data') sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) transformed_data = Array((sqrt_singular_values_inv * decomposition.Vt * data')') @@ -127,12 +130,12 @@ function GPObj( # Construct kernel: # Note that the kernels take the signal standard deviations on a # log scale as input. - rbf_len = log.(ones(size(inputs, 2))) + rbf_len = log.(ones(size(new_inputs, 2))) rbf_logstd = log(1.0) rbf = SEArd(rbf_len, rbf_logstd) kern=rbf - println("Using default squared exponential kernel:", kern) + println("Using default squared exponential kernel: ", kern) else kern=deepcopy(GPkernel) println("Using user-defined kernel",kern) @@ -162,23 +165,28 @@ function GPObj( GPkernel_i = deepcopy(kern) # inputs: N_samples x N_parameters # transformed_data: N_samples x N_data - + dat = dropdims(transformed_data[:, i]', dims=1) + #3dat = transformed_data[i,:] + # Zero mean function kmean = MeanZero() - m = GPE(inputs', - dropdims(transformed_data[:, i]', dims=1), + m = GPE(new_inputs', + dat, kmean, GPkernel_i, logstd_regularization_noise) - # we choose above to explicitly learn the WhiteKernel as opposed to using this - # in built functionality. + # we choose above to explicitly learn the WhiteKernel as opposed to using the + # in built noise=true functionality. + println("created GP", i) optimize!(m, noise=false) + println("optimized GP", i) push!(models, m) + print(m) end return GPObj{FT, typeof(package)}(inputs, - transformed_data, + data, input_mean, sqrt_inv_input_cov, models, @@ -216,7 +224,6 @@ function GPObj( # NB: svdfact() may be more efficient / provide positive definiteness information # stored as: svd.U * svd.S *svd.Vt decomposition = svd(truth_cov) - # Make sure inputs and data are arrays of type FT inputs = convert(Array{FT}, inputs) data = convert(Array{FT}, data) @@ -229,7 +236,9 @@ function GPObj( error("Cov not defined for 1d input; can't normalize") end sqrt_inv_input_cov = convert(Array{FT}, sqrt(inv(cov(inputs, dims=1)))) - inputs = (inputs .- input_mean) * sqrt_inv_input_cov + new_inputs = (inputs .- input_mean) * sqrt_inv_input_cov + else + new_inputs = inputs end # Transform the outputs @@ -243,7 +252,7 @@ function GPObj( const_value = 1.0 var_kern = ConstantKernel(constant_value=const_value, constant_value_bounds=(1e-05, 10000.0)) - rbf_len = ones(size(inputs, 2)) + rbf_len = ones(size(new_inputs, 2)) rbf = RBF(length_scale=rbf_len, length_scale_bounds=(1e-05, 10000.0)) @@ -266,7 +275,7 @@ function GPObj( #make the regularization small. We actually learn # total_noise = white_noise_level + regularization_noise magic_number = 1e-3 # magic_number << 1 - regularization_noise = 1e-3 + regularization_noise = 1e-3 else # when not learning noise, our SVD transformation implies the observational noise is the identity. regularization_noise = 1.0 @@ -275,18 +284,21 @@ function GPObj( for i in 1:size(data, 2) GPkernel_i = deepcopy(kern) - + # (previously) data= convert(Matrix{Float64},data') + # data=reshape(data[i,:],(size(data,2),1)) m = GaussianProcessRegressor(kernel=GPkernel_i, n_restarts_optimizer=10, alpha=regularization_noise) - - ScikitLearn.fit!(m, inputs, data[:, i]) + + println(m.kernel_) + ScikitLearn.fit!(m, new_inputs, transformed_data[:, i]) if i==1 println(m.kernel.hyperparameters) print("Completed training of: ") end print(i,", ") push!(models, m) + println(m.kernel_) end return GPObj{FT, typeof(package)}(inputs, data, @@ -320,7 +332,7 @@ function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, end M = length(gp.models) μσ2 = [predict_f(gp.models[i], new_inputs) for i in 1:M] - # Return mean(s) and standard deviation(s) + # Return mean(s) and variance(s) μ = reshape(vcat(first.(μσ2)...), size(new_inputs)[1], size(new_inputs)[2]) σ2 = reshape(vcat(last.(μσ2)...), @@ -354,11 +366,13 @@ function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, new_inputs = convert(Array{FT}, new_inputs') μσ2 = [predict_y(gp.models[i], new_inputs) for i in 1:M] - # Return mean(s) and standard deviation(s) + # Return mean(s) and variance(s) μ = reshape(vcat(first.(μσ2)...), size(new_inputs)[1], size(new_inputs)[2]) σ2 = reshape(vcat(last.(μσ2)...), size(new_inputs)[1], size(new_inputs)[2]) + + if transform_to_real # Revert back to the original (correlated) coordinate system # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP @@ -387,12 +401,13 @@ function predict(gp::GPObj{FT, SKLJL}, new_inputs::Array{FT}, transform_to_real) end M = length(gp.models) μσ = [gp.models[i].predict(new_inputs, return_std=true) for i in 1:M] - # Return mean(s) and standard deviation(s) - μ = reshape(vcat(first.(μσ2)...), + # Return mean(s) and standard deviations(s) + μ = reshape(vcat(first.(μσ)...), size(new_inputs)[1], size(new_inputs)[2]) - σ2 = reshape(vcat(last.(μσ2)...), + σ = reshape(vcat(last.(μσ)...), size(new_inputs)[1], size(new_inputs)[2]) - + σ2 = σ*σ + if transform_to_real # Revert back to the original (correlated) coordinate system # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP diff --git a/src/MCMC.jl b/src/MCMC.jl index e426aa10b..e8913772c 100644 --- a/src/MCMC.jl +++ b/src/MCMC.jl @@ -89,7 +89,7 @@ function MCMCObj( sqrt_singular_values_inv=Diagonal(1.0 ./ sqrt.(decomposition.S)) #diagonal matrix of 1/eigenvalues obs_sample=sqrt_singular_values_inv*decomposition.Vt*obs_sample else - println("Assuming independent outputs." + println("Assuming independent outputs.") end # first row is param_init From 2a272d78be2a70056301de51c32404ce4d821e2a Mon Sep 17 00:00:00 2001 From: odunbar Date: Wed, 8 Jul 2020 19:40:23 +0100 Subject: [PATCH 11/18] example working with single predictions! --- examples/GPEmulator/plot_GP.jl | 97 +++++++++++++++++++++++++--------- src/GPEmulator.jl | 34 ++++++++---- 2 files changed, 95 insertions(+), 36 deletions(-) diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl index 2302d0c7b..892b96f7e 100644 --- a/examples/GPEmulator/plot_GP.jl +++ b/examples/GPEmulator/plot_GP.jl @@ -36,18 +36,24 @@ Random.seed!(rng_seed) # The observables y are related to the parameters x by: # y = G(x1, x2) + η, # where G(x1, x2) := [sin(x1) + cos(x2), sin(x1) - cos(x2)], and η ~ N(0, Σ) -n = 100 # number of training points +n = 50 # number of training points p = 2 # input dim d = 2 # output dim -#X = 2.0 * π * rand(n, p) -x1 = range(0.0, stop=2*π, length=convert(Int64,sqrt(n))) -x2 = range(0.0, stop=2*π, length=convert(Int64,sqrt(n))) -X1, X2 = meshgrid(x1, x2) + +# Case 1: structured Grid +#x1 = range(0.0, stop=2*π, length=convert(Int64,sqrt(n))) +#x2 = range(0.0, stop=2*π, length=convert(Int64,sqrt(n))) +#X = hcat(X1[:], X2[:]) +#X1, X2 = meshgrid(x1, x2) + +# Case 2: random Grid +X = 2.0 * π * rand(n, p) + # Input for predict needs to be N_samples x N_parameters -X = hcat(X1[:], X2[:]) + μ = zeros(d) -Σ = 0.1 * [[0.8, 0.0] [0.0, 0.5]] # d x d +Σ = 0.1 * [[0.8, 0.2] [0.2, 0.5]] # d x d g1x = sin.(X[:, 1]) .+ cos.(X[:, 2]) g2x = sin.(X[:, 1]) .- cos.(X[:, 2]) @@ -55,51 +61,90 @@ noise_samples = rand(MvNormal(μ, Σ), n)' gx = [g1x g2x] Y = gx .+ noise_samples -truth_cov = Σ -sample_truth_cov = cov(noise_samples, dims=1) +# Case 1: input the true covariance +#input_cov = Σ +# Case 2: input a sample covariance approximation +input_cov = cov(noise_samples, dims=1) # Fit 2D Gaussian Process regression model # (To be precise, we fit two models, one that predicts y1 from x1 and x2, # and one that predicts y2 from x1 and x2) -gpobj = GPObj(X, Y, truth_cov, gppackage, GPkernel=nothing, - normalized=false, noise_learn=false, prediction_type=pred_type) +gpobj = GPObj(X, Y, input_cov, gppackage, GPkernel=nothing, + normalized=true, noise_learn=true, prediction_type=pred_type) # Plot mean and variance of the predicted observables y1 and y2 # For this, we generate test points on a x1-x2 grid -n_pts = 50 +n_pts = 200 x1 = range(0.0, stop=2*π, length=n_pts) x2 = range(0.0, stop=2*π, length=n_pts) X1, X2 = meshgrid(x1, x2) # Input for predict needs to be N_samples x N_parameters -inputs = hcat(X1[:], X2[:]) - -# Get predictions on the grid -gp_mean, gp_var = GPEmulator.predict(gpobj, inputs, transform_to_real=true) - +inputs = hcat(X1[:], X2[:]) + +gp_mean = zeros(size(inputs')) +gp_var = zeros(size(inputs')) +# Get predictions one-by-one on the grid +for pt in 1:size(inputs,1) + gp_mean_pt, gp_var_pt = GPEmulator.predict(gpobj, inputs[pt,:], transform_to_real=true) + #println(gp_mean_pt) + gp_mean[:,pt] = gp_mean_pt' + + gp_var[:,pt] = gp_var_pt' +end +# Make contour plots of the predictions for y_i in 1:d mean_grid = reshape(gp_mean[y_i, :], n_pts, n_pts) - p1 = plot(x1, x2, mean_grid, st=:surface, camera=(-30, 30), c=:cividis, - xlabel="x1", ylabel="x2", zlabel="mean of y"*string(y_i)) + p1 = plot(x1, x2, mean_grid, st=:contour, c=:cividis, + xlabel="x1", ylabel="x2") var_grid = reshape(gp_var[y_i, :], n_pts, n_pts) - p2 = plot(x1, x2, var_grid, st=:surface, camera=(-30, 30), c=:cividis, - xlabel="x1", ylabel="x2", zlabel="var of y"*string(y_i)) + p2 = plot(x1, x2, var_grid, st=:contour, c=:cividis, + xlabel="x1", ylabel="x2") plot(p1, p2, layout=(1, 2), legend=false) savefig("GP_test_y"*string(y_i)*".png") end -# Make plots of the true components of G(x1, x2) + +# Make contour plots of the true components of G(x1, x2) g1_true = sin.(inputs[:, 1]) .+ cos.(inputs[:, 2]) g1_grid = reshape(g1_true, n_pts, n_pts) -p3 = plot(x1, x2, g1_grid, st=:surface, camera=(-30, 30), c=:cividis, - xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") +p3 = plot(x1, x2, g1_grid, st=:contour, c=:cividis, + xlabel="x1", ylabel="x2") savefig("GP_test_true_g1.png") g2_true = sin.(inputs[:, 1]) .- cos.(inputs[:, 2]) g2_grid = reshape(g2_true, n_pts, n_pts) -p4 = plot(x1, x2, g2_grid, st=:surface, camera=(-30, 30), c=:cividis, - xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") +p4 = plot(x1, x2, g2_grid, st=:contout, c=:cividis, + xlabel="x1", ylabel="x2") savefig("GP_test_true_g2.png") + + + +# for y_i in 1:d +# mean_grid = reshape(gp_mean[y_i, :], n_pts, n_pts) +# p1 = plot(x1, x2, mean_grid, st=:surface, camera=(-30, 30), c=:cividis, +# xlabel="x1", ylabel="x2", zlabel="mean of y"*string(y_i)) + +# var_grid = reshape(gp_var[y_i, :], n_pts, n_pts) +# p2 = plot(x1, x2, var_grid, st=:surface, camera=(-30, 30), c=:cividis, +# xlabel="x1", ylabel="x2", zlabel="var of y"*string(y_i)) + +# plot(p1, p2, layout=(1, 2), legend=false) +# savefig("GP_test_y"*string(y_i)*".png") +# end + +# Make plots of the true components of G(x1, x2) +# g1_true = sin.(inputs[:, 1]) .+ cos.(inputs[:, 2]) +# g1_grid = reshape(g1_true, n_pts, n_pts) +# p3 = plot(x1, x2, g1_grid, st=:surface, camera=(-30, 30), c=:cividis, +# xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") +# savefig("GP_test_true_g1.png") + +# g2_true = sin.(inputs[:, 1]) .- cos.(inputs[:, 2]) +# g2_grid = reshape(g2_true, n_pts, n_pts) +# p4 = plot(x1, x2, g2_grid, st=:surface, camera=(-30, 30), c=:cividis, +# xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") +# savefig("GP_test_true_g2.png") diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 355ae0400..04270a7b7 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -18,6 +18,8 @@ export predict export GPJL, SKLJL export YType, FType +using Plots + """ GaussianProcessesPackage @@ -163,26 +165,35 @@ function GPObj( # Make a copy of GPkernel (because the kernel gets altered in # every iteration) GPkernel_i = deepcopy(kern) - # inputs: N_samples x N_parameters + # new_inputs: N_samples x N_parameters # transformed_data: N_samples x N_data - dat = dropdims(transformed_data[:, i]', dims=1) - #3dat = transformed_data[i,:] + GPinputs = new_inputs' + GPdata_i = dropdims(transformed_data[:, i]', dims=1) + + # for i=1:size(transformed_data,1) + # println(GPinputs[:,i]," ",GPdata_i[i]) + # end + + #dat = transformed_data[i,:] # Zero mean function kmean = MeanZero() - m = GPE(new_inputs', - dat, + m = GPE(GPinputs, + GPdata_i, kmean, GPkernel_i, logstd_regularization_noise) - + # we choose above to explicitly learn the WhiteKernel as opposed to using the # in built noise=true functionality. println("created GP", i) optimize!(m, noise=false) + # if i==1 + # plot(m, show = true) + # end println("optimized GP", i) push!(models, m) - print(m) + println(m.kernel) end return GPObj{FT, typeof(package)}(inputs, @@ -327,6 +338,7 @@ predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}; transform_to_real=false) whe function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, ::FType) where {FT} + if gp.normalized new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov end @@ -358,6 +370,9 @@ function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, end function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, ::YType) where {FT} + + #This just + new_inputs = convert(Array{FT}, new_inputs') if gp.normalized new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov end @@ -365,14 +380,13 @@ function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, # predicts columns of inputs so must be transposed new_inputs = convert(Array{FT}, new_inputs') μσ2 = [predict_y(gp.models[i], new_inputs) for i in 1:M] - + # Return mean(s) and variance(s) μ = reshape(vcat(first.(μσ2)...), - size(new_inputs)[1], size(new_inputs)[2]) + size(new_inputs)[1], size(new_inputs)[2]) σ2 = reshape(vcat(last.(μσ2)...), size(new_inputs)[1], size(new_inputs)[2]) - if transform_to_real # Revert back to the original (correlated) coordinate system # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP From 62838e3f69279c983b438db4f3f9ad46d64217ab Mon Sep 17 00:00:00 2001 From: Melanie Bieli Date: Thu, 9 Jul 2020 18:02:05 -0400 Subject: [PATCH 12/18] Get `plot_GP.jl` to work The example `plot_GP.jl` demonstrates how to plot the mean and variance of GP predictions (GP models are fitted to 2D input). The `predict` function in `GPEmulator.jl`has been modified such that it can predict multiple test points at once (passed as a matrix of size N_test_points x input_dim, i.e., each row corresponds to one point for which a prediction is made). Some dimensional checks have been added as well. --- Manifest.toml | 372 ++++++++++++++++++++++++--------- Project.toml | 1 + examples/GPEmulator/plot_GP.jl | 147 ++++++------- src/GPEmulator.jl | 295 +++++++++++++++----------- 4 files changed, 514 insertions(+), 301 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index e7eee69e8..6c9b19533 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -14,9 +14,9 @@ version = "0.3.3" [[Adapt]] deps = ["LinearAlgebra"] -git-tree-sha1 = "fd04049c7dd78cfef0b06cdc1f0f181467655712" +git-tree-sha1 = "0fac443759fa829ed8066db6cf1077d888bb6573" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "1.1.0" +version = "2.0.2" [[ArnoldiMethod]] deps = ["DelimitedFiles", "LinearAlgebra", "Random", "SparseArrays", "StaticArrays", "Test"] @@ -38,21 +38,21 @@ version = "3.5.0+3" [[ArrayInterface]] deps = ["LinearAlgebra", "Requires", "SparseArrays"] -git-tree-sha1 = "649c08a5a3a513f4662673d3777fe6ccb4df9f5d" +git-tree-sha1 = "0eccdcbe27fd6bd9cba3be31c67bdd435a21e865" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "2.8.7" +version = "2.9.1" [[ArrayLayouts]] deps = ["FillArrays", "LinearAlgebra"] -git-tree-sha1 = "89182776a99b69964e995cc2f1e37b5fc3476d56" +git-tree-sha1 = "a3254b3780a3544838ca0b7e23b1e9b06eb71bd8" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" -version = "0.3.4" +version = "0.3.5" [[BandedMatrices]] deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"] -git-tree-sha1 = "fd300e252fa1d96c75884cfa37fd6a5402c79d4b" +git-tree-sha1 = "eaf98fa821ab26c5825fa3a054db735755c335da" uuid = "aae01518-5342-5314-be14-df237901396f" -version = "0.15.12" +version = "0.15.15" [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" @@ -75,6 +75,12 @@ git-tree-sha1 = "3f2969de608af70db755cee9d4490a7294a6afc3" uuid = "764a87c0-6b3e-53db-9096-fe964310641d" version = "2.5.0" +[[Bzip2_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "3663bfffede2ef41358b6fc2e1d8a6d50b3c3904" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.6+2" + [[CEnum]] git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9" uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" @@ -92,17 +98,11 @@ git-tree-sha1 = "a6c17353ee38ddab30e73dcfaa1107752de724ec" uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597" version = "0.8.1" -[[ChainRules]] -deps = ["ChainRulesCore", "LinearAlgebra", "Reexport", "Requires", "Statistics"] -git-tree-sha1 = "85f130f2c5ce208a5a395b550802398d2fcc5ee6" -uuid = "082447d4-558c-5d27-93f4-14fc19e9eca2" -version = "0.6.4" - [[ChainRulesCore]] deps = ["MuladdMacro"] -git-tree-sha1 = "32e2c6e44d4fdd985b5688b5e85c1f6892cf3d15" +git-tree-sha1 = "4177411bef28d680948562abd25059dceb4d53ed" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "0.8.0" +version = "0.9.2" [[Cloudy]] deps = ["Coverage", "DifferentialEquations", "DocStringExtensions", "ForwardDiff", "HCubature", "LinearAlgebra", "Optim", "PyPlot", "SpecialFunctions", "TaylorSeries", "Test"] @@ -110,17 +110,23 @@ git-tree-sha1 = "7cb3ce76f5a63d17eb358b6030f7664931c75857" uuid = "9e3b23bb-e7cc-4b94-886c-65de2234ba87" version = "0.1.0" +[[ColorSchemes]] +deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"] +git-tree-sha1 = "7a15e3690529fd1042f0ab954dff7445b1efc8a5" +uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" +version = "3.9.0" + [[ColorTypes]] deps = ["FixedPointNumbers", "Random"] -git-tree-sha1 = "27eb374570946a02aa184ef5b403dabaa7380693" +git-tree-sha1 = "6e7aa35d0294f647bb9c985ccc34d4f5d371a533" uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" -version = "0.10.4" +version = "0.10.6" [[Colors]] deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Reexport"] -git-tree-sha1 = "1e9bba7984e78aa8cdeea7f9f7cc984ad4e4b1c7" +git-tree-sha1 = "5639e44833cfcf78c6a73fbceb4da75611d312cd" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" -version = "0.12.2" +version = "0.12.3" [[Combinatorics]] git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860" @@ -128,16 +134,16 @@ uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" version = "1.0.2" [[CommonSubexpressions]] -deps = ["Test"] -git-tree-sha1 = "efdaf19ab11c7889334ca247ff4c9f7c322817b0" +deps = ["MacroTools", "Test"] +git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" -version = "0.2.0" +version = "0.3.0" [[Compat]] deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "054993b6611376ddb40203e973e954fd9d1d1902" +git-tree-sha1 = "a6a8197ae253f2c1a22b2ae17c2dfaf5812c03aa" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.12.0" +version = "3.13.0" [[CompilerSupportLibraries_jll]] deps = ["Libdl", "Pkg"] @@ -162,6 +168,12 @@ git-tree-sha1 = "a2a6a5fea4d6f730ec4c18a76d27ec10e8ec1c50" uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" version = "1.0.0" +[[Contour]] +deps = ["StaticArrays"] +git-tree-sha1 = "81685fee51fc5168898e3cbd8b0f01506cd9148e" +uuid = "d38c429a-6771-53c6-b99e-75d170b6e991" +version = "0.5.4" + [[Coverage]] deps = ["CoverageTools", "HTTP", "JSON", "LibGit2", "MbedTLS"] git-tree-sha1 = "aad69c0c2d8080427f507c36b289c77fe465c9fe" @@ -187,15 +199,15 @@ version = "1.3.0" [[DataFrames]] deps = ["CategoricalArrays", "Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "Missings", "PooledArrays", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] -git-tree-sha1 = "02f08ae77249b7f6d4186b081a016fb7454c616f" +git-tree-sha1 = "d4436b646615928b634b37e99a3288588072f851" uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -version = "0.21.2" +version = "0.21.4" [[DataStructures]] deps = ["InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "be680f1ad03c0a03796aa3fda5a2180df7f83b46" +git-tree-sha1 = "edad9434967fdc0a2631a65d902228400642120c" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.17.18" +version = "0.17.19" [[DataValueInterfaces]] git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" @@ -218,9 +230,9 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" [[DiffEqBase]] deps = ["ArrayInterface", "ChainRulesCore", "ConsoleProgressMonitor", "DataStructures", "Distributed", "DocStringExtensions", "FunctionWrappers", "IterativeSolvers", "IteratorInterfaceExtensions", "LabelledArrays", "LinearAlgebra", "Logging", "LoggingExtras", "MuladdMacro", "Parameters", "Printf", "ProgressLogging", "RecipesBase", "RecursiveArrayTools", "RecursiveFactorization", "Requires", "Roots", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "TableTraits", "TerminalLoggers", "TreeViews", "ZygoteRules"] -git-tree-sha1 = "ae65fac7d9933f3d039c0296b5d41bf8c3d8f4ea" +git-tree-sha1 = "1c977b0c5219e1a4e3ae1ee16b2f29ccc3bd43b7" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.38.4" +version = "6.40.4" [[DiffEqCallbacks]] deps = ["DataStructures", "DiffEqBase", "ForwardDiff", "LinearAlgebra", "NLsolve", "OrdinaryDiffEq", "RecipesBase", "RecursiveArrayTools", "StaticArrays"] @@ -230,21 +242,21 @@ version = "2.13.3" [[DiffEqFinancial]] deps = ["DiffEqBase", "DiffEqNoiseProcess", "LinearAlgebra", "Markdown", "RandomNumbers"] -git-tree-sha1 = "f0c6f2b0b9fa463a90da06142e45ecf8e0b70bac" +git-tree-sha1 = "db08e0def560f204167c58fd0637298e13f58f73" uuid = "5a0ffddc-d203-54b0-88ba-2c03c0fc2e67" -version = "2.3.0" +version = "2.4.0" [[DiffEqJump]] deps = ["ArrayInterface", "Compat", "DataStructures", "DiffEqBase", "FunctionWrappers", "LinearAlgebra", "Parameters", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "StaticArrays", "Statistics", "TreeViews"] -git-tree-sha1 = "0dfb3aa84f687c78f697173665e304d28651dd55" +git-tree-sha1 = "c9baaba9b1ee1407473a2daac1ca0ffdcb72253d" uuid = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12" -version = "6.9.2" +version = "6.9.3" [[DiffEqNoiseProcess]] deps = ["DataStructures", "DiffEqBase", "Distributions", "LinearAlgebra", "PoissonRandom", "Random", "RandomNumbers", "RecipesBase", "RecursiveArrayTools", "Requires", "ResettableStacks", "StaticArrays", "Statistics"] -git-tree-sha1 = "fc9ba5c47246d1e6c15ae36ce9f5e67b6ffc06b7" +git-tree-sha1 = "474bba439ce886baab756744c54436d7628ef05e" uuid = "77a26b50-5914-5dd7-bc55-306e6241c503" -version = "4.2.0" +version = "4.3.0" [[DiffEqPhysics]] deps = ["DiffEqBase", "DiffEqCallbacks", "ForwardDiff", "LinearAlgebra", "Printf", "Random", "RecipesBase", "RecursiveArrayTools", "Reexport", "StaticArrays"] @@ -306,9 +318,9 @@ version = "0.24.11" [[ElasticArrays]] deps = ["Adapt"] -git-tree-sha1 = "d7d8924a01f5170d16cb24a24df20f027a45383a" +git-tree-sha1 = "b740f7f3fa87a2fa22c93582f2dbb9110c4ece6e" uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" -version = "1.2.2" +version = "1.2.3" [[ElasticPDMats]] deps = ["LinearAlgebra", "MacroTools", "PDMats", "Test"] @@ -317,10 +329,34 @@ uuid = "2904ab23-551e-5aed-883f-487f97af5226" version = "0.2.1" [[ExponentialUtilities]] -deps = ["LinearAlgebra", "Printf", "SparseArrays"] -git-tree-sha1 = "1672dedeacaab85345fd359ad56dde8fb5d48a45" +deps = ["LinearAlgebra", "Printf", "Requires", "SparseArrays"] +git-tree-sha1 = "91f7498b66205431fe3e35833cda97a22b1ab6a5" uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18" -version = "1.6.0" +version = "1.7.0" + +[[FFMPEG]] +deps = ["FFMPEG_jll"] +git-tree-sha1 = "c82bef6fc01e30d500f588cd01d29bdd44f1924e" +uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" +version = "0.3.0" + +[[FFMPEG_jll]] +deps = ["Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "LAME_jll", "LibVPX_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] +git-tree-sha1 = "0fa07f43e5609ea54848b82b4bb330b250e9645b" +uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" +version = "4.1.0+3" + +[[FFTW]] +deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"] +git-tree-sha1 = "14536c95939aadcee44014728a459d2fe3ca9acf" +uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +version = "1.2.2" + +[[FFTW_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "6c975cd606128d45d1df432fb812d6eb10fee00b" +uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" +version = "3.3.9+5" [[FastGaussQuadrature]] deps = ["LinearAlgebra", "SpecialFunctions"] @@ -330,15 +366,15 @@ version = "0.4.2" [[FillArrays]] deps = ["LinearAlgebra", "Random", "SparseArrays"] -git-tree-sha1 = "44f561e293987ffc84272cd3d2b14b0b93123d63" +git-tree-sha1 = "4783bbbeade37f2a8bd82af6c112510fde78e532" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "0.8.10" +version = "0.8.12" [[FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] -git-tree-sha1 = "fec7c2cb45c27071ef487fa7cae4fcac7509aa10" +git-tree-sha1 = "b02b6f6ea2c33f86a444f9cf132c1d1180a66cfd" uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" -version = "2.3.2" +version = "2.4.1" [[FixedPointNumbers]] git-tree-sha1 = "8fb797c37a3b7ced4327a05ac4ca0dd6a4f1ba92" @@ -353,9 +389,21 @@ version = "0.4.1" [[ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"] -git-tree-sha1 = "869540e4367122fbffaace383a5bdc34d6e5e5ac" +git-tree-sha1 = "1d090099fb82223abc48f7ce176d3f7696ede36d" uuid = "f6369f11-7733-5829-9624-2563aa707210" -version = "0.10.10" +version = "0.10.12" + +[[FreeType2_jll]] +deps = ["Bzip2_jll", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "7d900f32a3788d4eacac2bfa3bf5c770179c8afd" +uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7" +version = "2.10.1+2" + +[[FriBidi_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "2f56bee16bd0151de7b6a1eeea2ced190a2ad8d4" +uuid = "559328eb-81f9-559d-9380-de523a88c83c" +version = "1.0.5+3" [[FunctionWrappers]] git-tree-sha1 = "e4813d187be8c7b993cb7f85cbf2b7bfbaadc694" @@ -366,6 +414,12 @@ version = "1.1.1" deps = ["Random"] uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" +[[GR]] +deps = ["Base64", "DelimitedFiles", "HTTP", "JSON", "LinearAlgebra", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"] +git-tree-sha1 = "247adbd2b33c0c4b42efa20d1e807acf6312145f" +uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" +version = "0.50.1" + [[GaussianProcesses]] deps = ["Distances", "Distributions", "Documenter", "ElasticArrays", "ElasticPDMats", "FastGaussQuadrature", "ForwardDiff", "LinearAlgebra", "Optim", "PDMats", "Printf", "ProgressMeter", "Random", "RecipesBase", "ScikitLearnBase", "SpecialFunctions", "StaticArrays", "Statistics", "StatsFuns", "Zygote"] git-tree-sha1 = "9b33cd8be0fd2bc59c46911def57331a7d5f6e84" @@ -384,6 +438,18 @@ git-tree-sha1 = "62909c3eda8a25b5673a367d1ad2392ebb265211" uuid = "01680d73-4ee2-5a08-a1aa-533608c188bb" version = "0.3.0" +[[GeometryBasics]] +deps = ["IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"] +git-tree-sha1 = "119f32f9c2b497b49cd3f7f513b358b82660294c" +uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326" +version = "0.2.15" + +[[GeometryTypes]] +deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "StaticArrays"] +git-tree-sha1 = "34bfa994967e893ab2f17b864eec221b3521ba4d" +uuid = "4d00f742-c7ba-57c2-abde-4428a4b178cb" +version = "0.8.3" + [[HCubature]] deps = ["Combinatorics", "DataStructures", "LinearAlgebra", "QuadGK", "StaticArrays"] git-tree-sha1 = "fc4adde0a029bf69c0da3f7c25ca24b5e0e389ed" @@ -392,15 +458,15 @@ version = "1.4.0" [[HTTP]] deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"] -git-tree-sha1 = "ec87d5e2acbe1693789efbbe14f5ea7525758f71" +git-tree-sha1 = "eca61b35cdd8cd2fcc5eec1eda766424a995b02f" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "0.8.15" +version = "0.8.16" [[IRTools]] deps = ["InteractiveUtils", "MacroTools", "Test"] -git-tree-sha1 = "6875ae3cfcb9a50af80553d5cc825f406e8d13bc" +git-tree-sha1 = "90ee39f9beaaa186e4968417ea2b8ed5673c91c0" uuid = "7869d1d1-7146-5819-86e3-90919afe41df" -version = "0.4.0" +version = "0.3.3" [[Inflate]] git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c" @@ -413,6 +479,12 @@ git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8" uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f" version = "0.5.0" +[[IntelOpenMP_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "fb8e1c7a5594ba56f9011310790e03b5384998d6" +uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" +version = "2018.0.3+0" + [[InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" @@ -451,11 +523,17 @@ git-tree-sha1 = "8868479ff35ab98588ed0a529a9c2a4f8bda3bd6" uuid = "b14d175d-62b4-44ba-8fb7-3064adc8c3ec" version = "0.2.0" +[[LAME_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "221cc8998b9060677448cbb6375f00032554c4fd" +uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d" +version = "3.100.0+1" + [[LLVM]] deps = ["CEnum", "Libdl", "Printf", "Unicode"] -git-tree-sha1 = "72fc0a39d5899091ff2d4cdaa64cb5e4862cf813" +git-tree-sha1 = "a662366a5d485dee882077e8da3e1a95a86d097f" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "1.5.2" +version = "2.0.0" [[LaTeXStrings]] git-tree-sha1 = "de44b395389b84fd681394d4e8d39ef14e3a2ea8" @@ -464,9 +542,9 @@ version = "1.1.0" [[LabelledArrays]] deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"] -git-tree-sha1 = "f6def2c9c88908fdde81b9a39c236cf69de94450" +git-tree-sha1 = "5e04374019448f8509349948ab504f117e3b575a" uuid = "2ee39098-c373-598a-b85f-a56591580800" -version = "1.2.2" +version = "1.3.0" [[Latexify]] deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"] @@ -481,9 +559,14 @@ uuid = "1d6d02ad-be62-4b6b-8a6d-2f90e265016e" version = "0.1.2" [[LibGit2]] -deps = ["Printf"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" +[[LibVPX_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "e3549ca9bf35feb9d9d954f4c6a9032e92f46e7c" +uuid = "dd192d2f-8180-539f-9fb4-cc70b1dcf69a" +version = "1.8.1+1" + [[Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" @@ -513,9 +596,9 @@ version = "0.4.1" [[LoopVectorization]] deps = ["DocStringExtensions", "LinearAlgebra", "OffsetArrays", "SIMDPirates", "SLEEFPirates", "UnPack", "VectorizationBase"] -git-tree-sha1 = "59f7e9fddaae12967a0c0903aff2d06a8813e2b1" +git-tree-sha1 = "6f639e9cb9273c500f18779e4805646e6ac57d21" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.8.5" +version = "0.8.14" [[METIS_jll]] deps = ["Libdl", "Pkg"] @@ -523,6 +606,12 @@ git-tree-sha1 = "3f52ed323683398498ef163a45ce998f1ceca363" uuid = "d00139f3-1899-568f-a2f0-47f597d42d70" version = "5.1.0+4" +[[MKL_jll]] +deps = ["IntelOpenMP_jll", "Libdl", "Pkg"] +git-tree-sha1 = "0ce9a7fa68c70cf83c49d05d2c04d91b47404b08" +uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +version = "2020.1.216+0" + [[MLStyle]] git-tree-sha1 = "67f9a88611bc79f992aa705d9bbc833a2547dec7" uuid = "d8e11817-5142-5d16-987a-aa16d5891078" @@ -546,9 +635,14 @@ version = "1.0.2" [[MbedTLS_jll]] deps = ["Libdl", "Pkg"] -git-tree-sha1 = "c83f5a1d038f034ad0549f9ee4d5fac3fb429e33" +git-tree-sha1 = "a0cb0d489819fa7ea5f9fa84c7e7eba19d8073af" uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.16.0+2" +version = "2.16.6+1" + +[[Measures]] +git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f" +uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e" +version = "0.3.1" [[Missings]] deps = ["DataAPI"] @@ -606,15 +700,27 @@ uuid = "71a1bf82-56d0-4bbc-8a3c-48b961074391" version = "0.1.3" [[OffsetArrays]] -git-tree-sha1 = "ab697473e983a7499f463b696da8e8feb1b20ffd" +git-tree-sha1 = "4ba4cd84c88df8340da1c3e2d8dcb9d18dd1b53b" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.1.0" +version = "1.1.1" + +[[Ogg_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "59cf7a95bf5ac39feac80b796e0f39f9d69dc887" +uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051" +version = "1.3.4+0" [[OpenBLAS_jll]] deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"] -git-tree-sha1 = "1887096f6897306a4662f7c5af936da7d5d1a062" +git-tree-sha1 = "0c922fd9634e358622e333fc58de61f05a048492" uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.9+4" +version = "0.3.9+5" + +[[OpenSSL_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "7aaaded15bf393b5f34c2aad5b765c18d26cb495" +uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +version = "1.1.1+4" [[Optim]] deps = ["Compat", "FillArrays", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] @@ -622,10 +728,16 @@ git-tree-sha1 = "62054d469d3631960e3f472ceb8624be5b11c34d" uuid = "429524aa-4258-5aef-a3af-852621145aeb" version = "0.20.6" +[[Opus_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "002c18f222a542907e16c83c64a1338992da7e2c" +uuid = "91d4177d-7536-5919-b921-800302f37372" +version = "1.3.1+1" + [[OrderedCollections]] -git-tree-sha1 = "12ce190210d278e12644bcadf5b21cbdcf225cd3" +git-tree-sha1 = "293b70ac1780f9584c89268a6e2a560d938a7065" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.2.0" +version = "1.3.0" [[OrdinaryDiffEq]] deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "ExponentialUtilities", "FiniteDiff", "ForwardDiff", "GenericSVD", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "NLsolve", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"] @@ -653,14 +765,32 @@ version = "0.12.1" [[Parsers]] deps = ["Dates", "Test"] -git-tree-sha1 = "eb3e09940c0d7ae01b01d9291ebad7b081c844d3" +git-tree-sha1 = "10134f2ee0b1978ae7752c41306e131a684e1f06" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "1.0.5" +version = "1.0.7" [[Pkg]] -deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] +deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Test", "UUIDs"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +[[PlotThemes]] +deps = ["PlotUtils", "Requires", "Statistics"] +git-tree-sha1 = "c6f5ea535551b3b16835134697f0c65d06c94b91" +uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a" +version = "2.0.0" + +[[PlotUtils]] +deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"] +git-tree-sha1 = "e18e0e51ff07bf92bb7e06dcb9c082a4e125e20c" +uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043" +version = "1.0.5" + +[[Plots]] +deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "GeometryTypes", "JSON", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"] +git-tree-sha1 = "21d5dd05d230a98b97223d271561f1871989944c" +uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +version = "1.5.4" + [[PoissonRandom]] deps = ["Random", "Statistics", "Test"] git-tree-sha1 = "44d018211a56626288b5d3f8c6497d28c26dc850" @@ -714,9 +844,9 @@ version = "2.9.0" [[QuadGK]] deps = ["DataStructures", "LinearAlgebra"] -git-tree-sha1 = "dc84e810393cfc6294248c9032a9cdacc14a3db4" +git-tree-sha1 = "0ab8a09d4478ebeb99a706ecbf8634a65077ccdc" uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -version = "2.3.1" +version = "2.4.0" [[REPL]] deps = ["InteractiveUtils", "Markdown", "Sockets"] @@ -737,17 +867,23 @@ git-tree-sha1 = "54f8ceb165a0f6d083f0d12cb4996f5367c6edbc" uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" version = "1.0.1" +[[RecipesPipeline]] +deps = ["Dates", "PlotUtils", "RecipesBase"] +git-tree-sha1 = "d2a58b8291d1c0abae6a91489973f8a92bf5c04a" +uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c" +version = "0.1.11" + [[RecursiveArrayTools]] deps = ["ArrayInterface", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"] -git-tree-sha1 = "96e71928efa701fa5a6df0f88b51f05ceed70f2c" +git-tree-sha1 = "0ffe36b65f0fc4967a42a673c1a9ffa65724dee6" uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" -version = "2.4.4" +version = "2.5.0" [[RecursiveFactorization]] -deps = ["LinearAlgebra", "LoopVectorization"] -git-tree-sha1 = "09217cb106dd826de9960986207175b52e3035f2" +deps = ["LinearAlgebra", "LoopVectorization", "VectorizationBase"] +git-tree-sha1 = "4ca0bdad1d69abbd59c35af89a9a2ab6cd5ef0f1" uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" -version = "0.1.2" +version = "0.1.4" [[Reexport]] deps = ["Pkg"] @@ -781,24 +917,24 @@ version = "0.2.2+1" [[Roots]] deps = ["Printf"] -git-tree-sha1 = "c2f7348c55d1433d1cab0159b4d2c6d27af36fc4" +git-tree-sha1 = "7bd9d7eee602f0413f24c394386a18cb0d515f36" uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" -version = "1.0.2" +version = "1.0.3" [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" [[SIMDPirates]] deps = ["VectorizationBase"] -git-tree-sha1 = "74bf6ed250c21651955bdb36b2b12320374c49ae" +git-tree-sha1 = "dae629e96c1819d77882256e6cb29736f493bc30" uuid = "21efa798-c60a-11e8-04d3-e1a92915a26a" -version = "0.8.7" +version = "0.8.13" [[SLEEFPirates]] deps = ["Libdl", "SIMDPirates", "VectorizationBase"] -git-tree-sha1 = "3d97df9b38b3df1f118a203ac4a6c53b23265b4e" +git-tree-sha1 = "c750d618b7c8268a97e55c70e8c88e56080d30fa" uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa" -version = "0.5.1" +version = "0.5.4" [[SafeTestsets]] deps = ["Test"] @@ -825,6 +961,12 @@ uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" deps = ["Distributed", "Mmap", "Random", "Serialization"] uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" +[[Showoff]] +deps = ["Dates"] +git-tree-sha1 = "e032c9df551fb23c9f98ae1064de074111b7bc39" +uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f" +version = "0.3.1" + [[SimpleTraits]] deps = ["InteractiveUtils", "MacroTools"] git-tree-sha1 = "2ee666b24ab8be6a922f9d6c11a86e1a703a7dda" @@ -846,9 +988,9 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[SparseDiffTools]] deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "VertexSafeGraphs"] -git-tree-sha1 = "bfe68e0d914952932594b3c838f08463b0841037" +git-tree-sha1 = "93666e93899d995ec37abddde4811f533e49c074" uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804" -version = "1.8.0" +version = "1.9.1" [[SpecialFunctions]] deps = ["BinDeps", "BinaryProvider", "Libdl"] @@ -858,9 +1000,9 @@ version = "0.8.0" [[StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "5c06c0aeb81bef54aed4b3f446847905eb6cbda0" +git-tree-sha1 = "016d1e1a00fabc556473b07161da3d39726ded35" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "0.12.3" +version = "0.12.4" [[Statistics]] deps = ["LinearAlgebra", "SparseArrays"] @@ -890,6 +1032,12 @@ git-tree-sha1 = "2006e15cb3967afec91af1cb33b965eebea4f573" uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" version = "6.23.1" +[[StructArrays]] +deps = ["Adapt", "DataAPI", "Tables"] +git-tree-sha1 = "8099ed9fb90b6e754d6ba8c6ed8670f010eadca0" +uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" +version = "0.4.4" + [[SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" @@ -966,15 +1114,15 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" [[Unitful]] deps = ["ConstructionBase", "LinearAlgebra", "Random"] -git-tree-sha1 = "3714b55de06b11b2aa788b8643d6e91f13648be5" +git-tree-sha1 = "a061dada333813818aa7454f93c63a5cab6ea981" uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" -version = "1.2.1" +version = "1.3.0" [[VectorizationBase]] deps = ["CpuId", "LLVM", "Libdl", "LinearAlgebra"] -git-tree-sha1 = "bcadc352d9c81b0ef9ceebe822d30128b779f56b" +git-tree-sha1 = "bb72c58beab6c9e544851f5373fcd72f8f1f157a" uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" -version = "0.12.8" +version = "0.12.21" [[VersionParsing]] git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f" @@ -987,14 +1135,50 @@ git-tree-sha1 = "b9b450c99a3ca1cc1c6836f560d8d887bcbe356e" uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" version = "0.1.2" +[[Zlib_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "622d8b6dc0c7e8029f17127703de9819134d1b71" +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.11+14" + [[Zygote]] -deps = ["AbstractFFTs", "ArrayLayouts", "ChainRules", "FillArrays", "ForwardDiff", "Future", "IRTools", "InteractiveUtils", "LinearAlgebra", "MacroTools", "NNlib", "Random", "Requires", "Statistics", "ZygoteRules"] -git-tree-sha1 = "6d0f78976db6dbea9a36865efe068e6e2a5db6ed" +deps = ["DiffRules", "FFTW", "FillArrays", "ForwardDiff", "IRTools", "InteractiveUtils", "LinearAlgebra", "MacroTools", "NNlib", "NaNMath", "Random", "Requires", "SpecialFunctions", "Statistics", "ZygoteRules"] +git-tree-sha1 = "f8329b595c465caf3ca87c4f744e6041a4983e43" uuid = "e88e6eb3-aa80-5325-afca-941959d7151f" -version = "0.4.21" +version = "0.4.8" [[ZygoteRules]] deps = ["MacroTools"] git-tree-sha1 = "b3b4882cc9accf6731a08cc39543fbc6b669dca8" uuid = "700de1a5-db45-46bc-99cf-38207098b444" version = "0.2.0" + +[[libass_jll]] +deps = ["Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "027a304b2a90de84f690949a21f94e5ae0f92c73" +uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0" +version = "0.14.0+2" + +[[libfdk_aac_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "480c7ed04f68ea3edd4c757f5db5b6a0a4e0bd99" +uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280" +version = "0.1.6+2" + +[[libvorbis_jll]] +deps = ["Libdl", "Ogg_jll", "Pkg"] +git-tree-sha1 = "6a66f65b5275dfa799036c8a3a26616a0a271c4a" +uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a" +version = "1.3.6+4" + +[[x264_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "d89346fe63a6465a9f44e958ac0e3d366af90b74" +uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a" +version = "2019.5.25+2" + +[[x265_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "61324ad346b00a6e541896b94201c9426591e43a" +uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76" +version = "3.0.0+1" diff --git a/Project.toml b/Project.toml index 06c35bf68..2a5f22514 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324" diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl index 892b96f7e..e1344623a 100644 --- a/examples/GPEmulator/plot_GP.jl +++ b/examples/GPEmulator/plot_GP.jl @@ -3,14 +3,39 @@ using Random using GaussianProcesses using Distributions using Statistics -using Plots; pyplot(size = (1000, 400)) +using Plots; pyplot(size = (1500, 700)) +Plots.scalefontsizes(1.3) using LinearAlgebra using CalibrateEmulateSample.GPEmulator ############################################################################### # # -# This examples shows how to fit a Gaussian Process regression model to # -# 2D input data, and how to plot the mean and variance # +# This examples shows how to fit a Gaussian Process regression model # +# using GPEmulator, and how to plot the mean and variance # +# # +# Training points: {(x_i, y_i)} (i=1,...,n), where x_i ∈ ℝ ² and y_i ∈ ℝ ² # +# The y_i are assumed to be related to the x_i by: # +# y_i = G(x_i) + η, where G is a map from ℝ ² to ℝ ², and η is noise drawn # +# from a 2D normal distribution with known mean μ and covariance Σ # +# # +# Two Gaussian Process models are fit, one that predicts y_i[1] from x_i # +# and one that predicts y_i[2] from x_i. Fitting separate models for # +# y_i[1] and y_i[2] requires the outputs to be independent - since this # +# cannot be assumed to be the case a priori, the training data are # +# decorrelated by perfoming an SVD decomposition on the noise # +# covariance Σ, and each model is then trained in the decorrelated space. # +# # +# The decorrelation is done automatically when instantiating a `GPObj`, # +# but the user can choose to have the `predict` function return its # +# predictions in the original space (by setting transform_to_real=true) # +# # +# The GPEmulator module can be used as a standalone module to fit # +# Gaussian Process regression models, but it was originally designed as # +# the "Emulate" step of the "Calibrate-Emulate-Sample" framework # +# developed at CliMA. # +# Further reading on Calibrate-Emulate-Sample: # +# Cleary et al., 2020 # +# https://arxiv.org/abs/2001.03689 # # # ############################################################################### @@ -22,13 +47,11 @@ function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T return gx, gy end -gppackage = GPJL() -pred_type = YType() - -# Seed for pseudo-random number generator rng_seed = 41 Random.seed!(rng_seed) +gppackage = GPJL() +pred_type = YType() # Generate training data (x-y pairs, where x ∈ ℝ ᵖ, y ∈ ℝ ᵈ) # x = [x1, x2]: inputs/predictors/features/parameters @@ -36,40 +59,33 @@ Random.seed!(rng_seed) # The observables y are related to the parameters x by: # y = G(x1, x2) + η, # where G(x1, x2) := [sin(x1) + cos(x2), sin(x1) - cos(x2)], and η ~ N(0, Σ) -n = 50 # number of training points +n = 50 # number of training points p = 2 # input dim d = 2 # output dim -# Case 1: structured Grid -#x1 = range(0.0, stop=2*π, length=convert(Int64,sqrt(n))) -#x2 = range(0.0, stop=2*π, length=convert(Int64,sqrt(n))) -#X = hcat(X1[:], X2[:]) -#X1, X2 = meshgrid(x1, x2) - -# Case 2: random Grid X = 2.0 * π * rand(n, p) +# G(x1, x2) +g1x = sin.(X[:, 1]) .+ cos.(X[:, 2]) +g2x = sin.(X[:, 1]) .- cos.(X[:, 2]) +gx = [g1x g2x] -# Input for predict needs to be N_samples x N_parameters - +# Add noise η μ = zeros(d) Σ = 0.1 * [[0.8, 0.2] [0.2, 0.5]] # d x d - -g1x = sin.(X[:, 1]) .+ cos.(X[:, 2]) -g2x = sin.(X[:, 1]) .- cos.(X[:, 2]) noise_samples = rand(MvNormal(μ, Σ), n)' -gx = [g1x g2x] -Y = gx .+ noise_samples -# Case 1: input the true covariance -#input_cov = Σ -# Case 2: input a sample covariance approximation +# y = G(x) + η +Y = gx .+ noise_samples input_cov = cov(noise_samples, dims=1) - # Fit 2D Gaussian Process regression model -# (To be precise, we fit two models, one that predicts y1 from x1 and x2, +# (To be precise: We fit two models, one that predicts y1 from x1 and x2, # and one that predicts y2 from x1 and x2) +# Setting GPkernel=nothing leads to the creation of a default squared +# exponential kernel. +# Setting noise_learn=true leads to the addition of white noise to the +# kernel gpobj = GPObj(X, Y, input_cov, gppackage, GPkernel=nothing, normalized=true, noise_learn=true, prediction_type=pred_type) @@ -79,72 +95,39 @@ n_pts = 200 x1 = range(0.0, stop=2*π, length=n_pts) x2 = range(0.0, stop=2*π, length=n_pts) X1, X2 = meshgrid(x1, x2) - -# Input for predict needs to be N_samples x N_parameters +# Input for predict has to be of size N_samples x input_dim inputs = hcat(X1[:], X2[:]) -gp_mean = zeros(size(inputs')) -gp_var = zeros(size(inputs')) -# Get predictions one-by-one on the grid -for pt in 1:size(inputs,1) - gp_mean_pt, gp_var_pt = GPEmulator.predict(gpobj, inputs[pt,:], transform_to_real=true) - #println(gp_mean_pt) - gp_mean[:,pt] = gp_mean_pt' - - gp_var[:,pt] = gp_var_pt' -end -# Make contour plots of the predictions for y_i in 1:d - mean_grid = reshape(gp_mean[y_i, :], n_pts, n_pts) - p1 = plot(x1, x2, mean_grid, st=:contour, c=:cividis, - xlabel="x1", ylabel="x2") - - var_grid = reshape(gp_var[y_i, :], n_pts, n_pts) - p2 = plot(x1, x2, var_grid, st=:contour, c=:cividis, - xlabel="x1", ylabel="x2") + gp_mean, gp_var = GPEmulator.predict(gpobj, + inputs, + transform_to_real=true) + mean_grid = reshape(gp_mean[:, y_i], n_pts, n_pts) + p1 = plot(x1, x2, mean_grid, st=:surface, camera=(-30, 30), c=:cividis, + xlabel="x1", ylabel="x2", zlabel="mean of y"*string(y_i), + zguidefontrotation=90) + + var_grid = reshape(gp_var[:, y_i], n_pts, n_pts) + p2 = plot(x1, x2, var_grid, st=:surface, camera=(-30, 30), c=:cividis, + xlabel="x1", ylabel="x2", zlabel="var of y"*string(y_i), + zguidefontrotation=90) plot(p1, p2, layout=(1, 2), legend=false) - savefig("GP_test_y"*string(y_i)*".png") + savefig("GP_test_y"*string(y_i)*"_predictions.png") end - -# Make contour plots of the true components of G(x1, x2) +# Make plots of the true components of G(x1, x2) g1_true = sin.(inputs[:, 1]) .+ cos.(inputs[:, 2]) g1_grid = reshape(g1_true, n_pts, n_pts) -p3 = plot(x1, x2, g1_grid, st=:contour, c=:cividis, - xlabel="x1", ylabel="x2") +p3 = plot(x1, x2, g1_grid, st=:surface, camera=(-30, 30), c=:cividis, + xlabel="x1", ylabel="x2", zlabel="sin(x1) + cos(x2)", + zguidefontrotation=90) savefig("GP_test_true_g1.png") g2_true = sin.(inputs[:, 1]) .- cos.(inputs[:, 2]) g2_grid = reshape(g2_true, n_pts, n_pts) -p4 = plot(x1, x2, g2_grid, st=:contout, c=:cividis, - xlabel="x1", ylabel="x2") +p4 = plot(x1, x2, g2_grid, st=:surface, camera=(-30, 30), c=:cividis, + xlabel="x1", ylabel="x2", zlabel="sin(x1) - cos(x2)", + zguidefontrotation=90) savefig("GP_test_true_g2.png") - - - -# for y_i in 1:d -# mean_grid = reshape(gp_mean[y_i, :], n_pts, n_pts) -# p1 = plot(x1, x2, mean_grid, st=:surface, camera=(-30, 30), c=:cividis, -# xlabel="x1", ylabel="x2", zlabel="mean of y"*string(y_i)) - -# var_grid = reshape(gp_var[y_i, :], n_pts, n_pts) -# p2 = plot(x1, x2, var_grid, st=:surface, camera=(-30, 30), c=:cividis, -# xlabel="x1", ylabel="x2", zlabel="var of y"*string(y_i)) - -# plot(p1, p2, layout=(1, 2), legend=false) -# savefig("GP_test_y"*string(y_i)*".png") -# end - -# Make plots of the true components of G(x1, x2) -# g1_true = sin.(inputs[:, 1]) .+ cos.(inputs[:, 2]) -# g1_grid = reshape(g1_true, n_pts, n_pts) -# p3 = plot(x1, x2, g1_grid, st=:surface, camera=(-30, 30), c=:cividis, -# xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") -# savefig("GP_test_true_g1.png") - -# g2_true = sin.(inputs[:, 1]) .- cos.(inputs[:, 2]) -# g2_grid = reshape(g2_true, n_pts, n_pts) -# p4 = plot(x1, x2, g2_grid, st=:surface, camera=(-30, 30), c=:cividis, -# xlabel="x1", ylabel="x2", zlabel="sin(x1) + 0.5*x2") -# savefig("GP_test_true_g2.png") +Plots.scalefontsizes(1/1.3) diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 04270a7b7..37ffd7983 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -54,13 +54,13 @@ $(DocStringExtensions.FIELDS) """ struct GPObj{FT<:AbstractFloat, GPM} - "training inputs/parameters; N_samples x N_parameters" + "training inputs/parameters; N_samples x input_dim" inputs::Array{FT} - "training data/targets; N_samples x N_data" + "training data/targets; N_samples x output_dim" data::Array{FT} - "mean of input; 1 x N_parameters" + "mean of input; 1 x input_dim" input_mean::Array{FT} - "square root of the inverse of the input covariance matrix; N_parameters x N_parameters" + "square root of the inverse of the input covariance matrix; input_dim x input_dim" sqrt_inv_input_cov::Union{Nothing, Array{FT, 2}} "the Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs" models::Vector @@ -78,7 +78,7 @@ end """ GPObj(inputs, data, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, svdflag=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} -Inputs and data of size N_samples x N_parameters (both arrays will be transposed in the construction of the GPObj) +Inputs and data of size N_samples x input_dim (both arrays will be transposed in the construction of the GPObj) - `GPkernel` - GaussianProcesses kernel object. If not supplied, a default kernel is used. The default kernel is the sum of a Squared @@ -95,102 +95,103 @@ function GPObj( noise_learn::Bool=true, prediction_type::PredictionType=YType()) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject} - - #FT = eltype(data) - models = Any[] + # Consistency check + err = "Number of inputs is not equal to the number of data." + size(inputs, 1) == size(data, 1) || throw(ArgumentError(err)) - # create an SVD decomposition of the covariance: - # NB: svdfact() may be more efficient / provide positive definiteness information - # stored as: svd.U * svd.S *svd.Vt - decomposition = svd(truth_cov) + # Initialize vector of GP models + models = Any[] + # Number of models (We are fitting one model per output dimension) + N_models = size(data, 2) # Make sure inputs and data are arrays of type FT inputs = convert(Array{FT}, inputs) data = convert(Array{FT}, data) - # Transform the inputs + # Normalize the inputs if normalized==true input_mean = reshape(mean(inputs, dims=1), 1, :) sqrt_inv_input_cov = nothing if normalized + # Normalize and transpose (since the inputs have to be of size + # input_dim x N_samples to pass to GPE()) if ndims(inputs) == 1 error("Cov not defined for 1d input; can't normalize") end sqrt_inv_input_cov = convert(Array{FT}, sqrt(inv(cov(inputs, dims=1)))) - new_inputs = (inputs .- input_mean) * sqrt_inv_input_cov + GPinputs = convert(Array, ((inputs.-input_mean) * sqrt_inv_input_cov)') else - new_inputs = inputs + # Only transpose + GPinputs = convert(Array, inputs') end - # Transform the outputs - #data = convert(Matrix{Float64},data') + # Transform the outputs (transformed_data will also come out transposed from + # this transformation, i.e., with size output_dim x N_samples. This is also + # the size that will be required when passing the data to GPE(), so we leave + # transformed_data in that size. + # Create an SVD decomposition of the covariance: + # NB: svdfact() may be more efficient / provide positive definiteness information + # stored as: svd.U * svd.S *svd.Vt + decomposition = svd(truth_cov) sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - transformed_data = Array((sqrt_singular_values_inv * decomposition.Vt * data')') + transformed_data = sqrt_singular_values_inv * decomposition.Vt * data' # Use a default kernel unless a kernel was supplied to GPObj if GPkernel==nothing - println("Using default squared exponential kernel, learning lengthscale and variance parameters") + println("Using default squared exponential kernel, learning length scale and variance parameters") # Construct kernel: # Note that the kernels take the signal standard deviations on a # log scale as input. - rbf_len = log.(ones(size(new_inputs, 2))) + rbf_len = log.(ones(size(GPinputs, 1))) rbf_logstd = log(1.0) rbf = SEArd(rbf_len, rbf_logstd) - - kern=rbf + kern = rbf println("Using default squared exponential kernel: ", kern) else - kern=deepcopy(GPkernel) - println("Using user-defined kernel",kern) + kern = deepcopy(GPkernel) + println("Using user-defined kernel", kern) end if noise_learn - # regularize with white noise + # Add white noise to kernel white_logstd = log(1.0) white = Noise(white_logstd) - - # add to kernel construct kernel kern = kern + white println("Learning additive white noise") - #make the regularization small. We actually learn - # total_noise = white_logstd + logstd_regularization_noise) + # Make the regularization small. We actually learn + # total_noise = white_logstd + logstd_regularization_noise magic_number = 1e-3 # magic_number << 1 - logstd_regularization_noise = log(sqrt(1e-3)) + logstd_regularization_noise = log(sqrt(magic_number)) + else - # when not learning noise, our SVD transformation implies the observational noise is the identity. + # When not learning noise, our SVD transformation implies the + # observational noise is the identity. logstd_regularization_noise = log(sqrt(1.0)) end - for i in 1:size(transformed_data, 2) - # Make a copy of GPkernel (because the kernel gets altered in - # every iteration) + for i in 1:N_models + # Make a copy of the kernel (because it gets altered in every + # iteration) GPkernel_i = deepcopy(kern) - # new_inputs: N_samples x N_parameters - # transformed_data: N_samples x N_data - GPinputs = new_inputs' - GPdata_i = dropdims(transformed_data[:, i]', dims=1) - - # for i=1:size(transformed_data,1) - # println(GPinputs[:,i]," ",GPdata_i[i]) - # end - - #dat = transformed_data[i,:] + GPdata_i = transformed_data[i, :] + # GPE() arguments: + # GPinputs: (input_dim x N_samples) + # GPdata_i: (N_samples,) # Zero mean function kmean = MeanZero() + + # Instantiate GP model m = GPE(GPinputs, GPdata_i, kmean, GPkernel_i, logstd_regularization_noise) - # we choose above to explicitly learn the WhiteKernel as opposed to using the - # in built noise=true functionality. + # We choose above to explicitly learn the WhiteKernel as opposed to + # using the in built noise=true functionality. println("created GP", i) optimize!(m, noise=false) - # if i==1 - # plot(m, show = true) - # end println("optimized GP", i) push!(models, m) println(m.kernel) @@ -208,11 +209,10 @@ function GPObj( end - """ GPObj(inputs, data, package::SKLJL, GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true) where {K<:Kernel, KPy<:PyObject} -Inputs and data of size N_samples x N_parameters (both arrays will be transposed in the construction of the GPObj) +Inputs and data of size N_samples x input_dim (both arrays will be transposed in the construction of the GPObj) - `GPkernel` - GaussianProcesses or ScikitLearn kernel object. If not supplied, a default kernel is used. The default kernel is the sum of a @@ -228,18 +228,23 @@ function GPObj( noise_learn::Bool=true, prediction_type::PredictionType=YType()) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject} - #FT = eltype(data) - models = Any[] + # Consistency check + err = "Number of inputs is not equal to the number of data." + size(inputs, 1) == size(data, 1) || throw(ArgumentError(err)) - # create an SVD decomposition of the covariance: - # NB: svdfact() may be more efficient / provide positive definiteness information - # stored as: svd.U * svd.S *svd.Vt - decomposition = svd(truth_cov) # Make sure inputs and data are arrays of type FT inputs = convert(Array{FT}, inputs) data = convert(Array{FT}, data) - # Transform the inputs + # Initialize vector of GP models + models = Any[] + # Number of models (We are fitting one model per output dimension) + N_models = size(data, 2) + + # Normalize the inputs if normalized==true + # Note that contrary to the GaussianProcesses.jl (GPJL) GPE, the + # ScikitLearn (SKLJL) GaussianProcessRegressor requires inputs to be of + # size N_samples x input_dim, so no transposition is needed here input_mean = reshape(mean(inputs, dims=1), 1, :) sqrt_inv_input_cov = nothing if normalized @@ -247,26 +252,27 @@ function GPObj( error("Cov not defined for 1d input; can't normalize") end sqrt_inv_input_cov = convert(Array{FT}, sqrt(inv(cov(inputs, dims=1)))) - new_inputs = (inputs .- input_mean) * sqrt_inv_input_cov + GPinputs = (inputs .- input_mean) * sqrt_inv_input_cov else - new_inputs = inputs + GPinputs = inputs end # Transform the outputs + # Create an SVD decomposition of the covariance: + # NB: svdfact() may be more efficient / provide positive definiteness information + # stored as: svd.U * svd.S *svd.Vt + decomposition = svd(truth_cov) sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - #transformed_data = sqrt_singular_values_inv * decomposition.Vt * data - transformed_data = Array((sqrt_singular_values_inv * decomposition.Vt * data')') + transformed_data = sqrt_singular_values_inv * decomposition.Vt * data' if GPkernel==nothing - println("Using default squared exponential kernel, learning lengthscale and variance parameters") - #Create default squared exponential kernel + println("Using default squared exponential kernel, learning length scale and variance parameters") + # Create default squared exponential kernel const_value = 1.0 var_kern = ConstantKernel(constant_value=const_value, - constant_value_bounds=(1e-05, 10000.0)) - rbf_len = ones(size(new_inputs, 2)) + constant_value_bounds=(1e-05, 10000.0)) + rbf_len = ones(size(GPinputs, 2)) rbf = RBF(length_scale=rbf_len, length_scale_bounds=(1e-05, 10000.0)) - - kern = var_kern * rbf println("Using default squared exponential kernel:", kern) else @@ -275,34 +281,35 @@ function GPObj( end if noise_learn - # regularize with white noise + # Add white noise to kernel white_noise_level = 1.0 white = WhiteKernel(noise_level=white_noise_level, noise_level_bounds=(1e-05, 10.0)) - kern = kern + white println("Learning additive white noise") - #make the regularization small. We actually learn + # Make the regularization small. We actually learn # total_noise = white_noise_level + regularization_noise magic_number = 1e-3 # magic_number << 1 regularization_noise = 1e-3 else - # when not learning noise, our SVD transformation implies the observational noise is the identity. + # When not learning noise, our SVD transformation implies the + # observational noise is the identity. regularization_noise = 1.0 end - - for i in 1:size(data, 2) + for i in 1:N_models GPkernel_i = deepcopy(kern) - # (previously) data= convert(Matrix{Float64},data') - # data=reshape(data[i,:],(size(data,2),1)) + GPdata_i = transformed_data[i, :] m = GaussianProcessRegressor(kernel=GPkernel_i, n_restarts_optimizer=10, alpha=regularization_noise) println(m.kernel_) - ScikitLearn.fit!(m, new_inputs, transformed_data[:, i]) + # ScikitLearn.fit! arguments: + # GPinputs: (N_samples x input_dim) + # GPdata_i: (N_samples,) + ScikitLearn.fit!(m, GPinputs, GPdata_i) if i==1 println(m.kernel.hyperparameters) print("Completed training of: ") @@ -325,121 +332,159 @@ end """ - predict(gp::GPObj{FT,GPJL}, new_inputs::Array{FT}) where {FT} = predict(gp, new_inputs, gp.prediction_type) + predict(gp::GPObj{FT,GPJL}, new_inputs::Array{FT, 2}) where {FT} = predict(gp, new_inputs, gp.prediction_type) Evaluate the GP model(s) at new inputs. - `gp` - a GPObj - - `new_inputs` - inputs for which GP model(s) is/are evaluated; N_new_inputs x N_parameters + - `new_inputs` - inputs for which GP model(s) is/are evaluated; N_new_inputs x input_dim Note: If gp.normalized == true, the new inputs are normalized prior to the prediction """ -predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}; transform_to_real=false) where {FT} = predict(gp, new_inputs, transform_to_real, gp.prediction_type) - -function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, - ::FType) where {FT} +predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}; transform_to_real=false) where {FT} = predict(gp, new_inputs, transform_to_real, gp.prediction_type) + +function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}, transform_to_real, ::FType) where {FT} + # Check if the size of new_inputs is consistent with the GP model's input + # dimension. + # new_inputs should be of size N_new_inputs x input_dim, so input_dim + # has to equal the input dimension of the GP model(s). + input_dim = size(gp.inputs, 2) + output_dim = size(gp.data, 2) + err = "GP object and input observations do not have consistent dimensions" + size(new_inputs, 2) == input_dim || throw(ArgumentError(err)) - if gp.normalized new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov end + M = length(gp.models) μσ2 = [predict_f(gp.models[i], new_inputs) for i in 1:M] # Return mean(s) and variance(s) - μ = reshape(vcat(first.(μσ2)...), - size(new_inputs)[1], size(new_inputs)[2]) - σ2 = reshape(vcat(last.(μσ2)...), - size(new_inputs)[1], size(new_inputs)[2]) + μ = reshape(vcat(first.(μσ2)...), size(new_inputs, 1), output_dim) + σ2 = reshape(vcat(last.(μσ2)...), size(new_inputs, 1), output_dim) + if transform_to_real # Revert back to the original (correlated) coordinate system - # We created meanvGP = Dinv*Vt*meanv, so meanv = V*D*meanvGP - μ = gp.decomposition.V * gp.sqrt_singular_values * μ - temp = zeros(size(new_inputs)) + # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP + transformed_μT= gp.decomposition.V * gp.sqrt_singular_values * μ' + # Get the means back into size N_new_inputs x output_dim + transformed_μ = convert(Array, transformed_μT') + transformed_σ2 = zeros(size(σ2)) # Back transformation of covariance - for j in 1:size(new_inputs, 2) - σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[:,j]) * gp.sqrt_singular_values * gp.decomposition.Vt + for j in 1:size(σ2, 1) + σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[j, :]) * gp.sqrt_singular_values * gp.decomposition.Vt # Extract diagonal elements from σ2_j and add them to σ2 as column # vectors # TODO: Should we have an optional input argument `full_cov` - if # set to true, the full covariance matrix is returned instead of # only variances (default: false)? - temp[:, j] = diag(σ2_j) + transformed_σ2[j, :] = diag(σ2_j) end - σ2 = temp + return transformed_μ, transformed_σ2 + + else + # No back transformation needed + return μ, σ2 end - return μ, σ2 end -function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT}, transform_to_real, ::YType) where {FT} +function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}, transform_to_real, ::YType) where {FT} + # Check if the size of new_inputs is consistent with the GP model's input + # dimension. + # new_inputs should be of size N_new_inputs x input_dim, so input_dim + # has to equal the input dimension of the GP model(s). + input_dim = size(gp.inputs, 2) + output_dim = size(gp.data, 2) + err = "GP object and input observations do not have consistent dimensions" + size(new_inputs, 2) == input_dim || throw(ArgumentError(err)) - #This just - new_inputs = convert(Array{FT}, new_inputs') if gp.normalized new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov end + M = length(gp.models) - # predicts columns of inputs so must be transposed - new_inputs = convert(Array{FT}, new_inputs') - μσ2 = [predict_y(gp.models[i], new_inputs) for i in 1:M] + # Predicts columns of inputs so must be transposed to size + # input_dim x N_new_inputs + new_inputs_transp = convert(Array{FT}, new_inputs') + μσ2 = [predict_y(gp.models[i], new_inputs_transp) for i in 1:M] # Return mean(s) and variance(s) - μ = reshape(vcat(first.(μσ2)...), - size(new_inputs)[1], size(new_inputs)[2]) - σ2 = reshape(vcat(last.(μσ2)...), - size(new_inputs)[1], size(new_inputs)[2]) + μ = reshape(vcat(first.(μσ2)...), size(new_inputs, 1), output_dim) + σ2 = reshape(vcat(last.(μσ2)...), size(new_inputs, 1), output_dim) if transform_to_real # Revert back to the original (correlated) coordinate system # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP - μ = gp.decomposition.V * gp.sqrt_singular_values * μ - temp = zeros(size(σ2)) + transformed_μT= gp.decomposition.V * gp.sqrt_singular_values * μ' + # Get the means back into size N_new_inputs x output_dim + transformed_μ = convert(Array, transformed_μT') + transformed_σ2 = zeros(size(σ2)) # Back transformation of covariance - for j in 1:size(σ2, 2) - σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[:, j]) * gp.sqrt_singular_values * gp.decomposition.Vt + for j in 1:size(σ2, 1) + σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[j, :]) * gp.sqrt_singular_values * gp.decomposition.Vt # Extract diagonal elements from σ2_j and add them to σ2 as column # vectors # TODO: Should we have an optional input argument `full_cov` - if # set to true, the full covariance matrix is returned instead of # only variances (default: false)? - temp[:, j] = diag(σ2_j) + transformed_σ2[j, :] = diag(σ2_j) end - σ2 = temp + return transformed_μ, transformed_σ2 + + else + # No back transformation needed + return μ, σ2 end - return μ, σ2 end -function predict(gp::GPObj{FT, SKLJL}, new_inputs::Array{FT}, transform_to_real) where {FT} +function predict(gp::GPObj{FT, SKLJL}, new_inputs::Array{FT, 2}, transform_to_real) where {FT} + + # Check if the size of new_inputs is consistent with the GP model's input + # dimension. + # new_inputs should be of size N_new_inputs x input_dim, so input_dim + # has to equal the input dimension of the GP model(s). + input_dim = size(gp.inputs, 2) + output_dim = size(gp.data, 2) + err = "GP object and input observations do not have consistent dimensions" + size(new_inputs, 2) == input_dim || throw(ArgumentError(err)) if gp.normalized new_inputs = (new_inputs .- gp.input_mean) * gp.sqrt_inv_input_cov end + M = length(gp.models) + # Predicts rows of inputs; no need to transpose μσ = [gp.models[i].predict(new_inputs, return_std=true) for i in 1:M] # Return mean(s) and standard deviations(s) - μ = reshape(vcat(first.(μσ)...), - size(new_inputs)[1], size(new_inputs)[2]) - σ = reshape(vcat(last.(μσ)...), - size(new_inputs)[1], size(new_inputs)[2]) + μ = reshape(vcat(first.(μσ)...), size(new_inputs)[1], output_dim) + σ = reshape(vcat(last.(μσ)...), size(new_inputs)[1], output_dim) σ2 = σ*σ if transform_to_real # Revert back to the original (correlated) coordinate system # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP - μ = gp.decomposition.V * gp.sqrt_singular_values * μ - temp = zeros(size(σ2)) + transformed_μT= gp.decomposition.V * gp.sqrt_singular_values * μ' + # Get the means back into size N_new_inputs x output_dim + transformed_μ = convert(Array, transformed_μT') + println("size transformed μ") + println(size(transformed_μ)) + transformed_σ2 = zeros(size(σ2)) # Back transformation of covariance - for j in 1:size(σ2, 2) - σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[:, j]) * gp.sqrt_singular_values * gp.decomposition.Vt + for j in 1:size(σ2, 1) + σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[j, :]) * gp.sqrt_singular_values * gp.decomposition.Vt # Extract diagonal elements from σ2_j and add them to σ2 as column # vectors # TODO: Should we have an optional input argument `full_cov` - if # set to true, the full covariance matrix is returned instead of # only variances (default: false)? - temp[:, j] = diag(σ2_j) + transformed_σ2[j, :] = diag(σ2_j) end - σ2 = temp + return transformed_μ, transformed_σ2 + + else + # No back transformation needed + return μ, σ2 end - return μ, σ2 end From 4a61651713836dc7f07b2542f3a465635704ae7a Mon Sep 17 00:00:00 2001 From: Melanie Bieli Date: Tue, 14 Jul 2020 01:03:24 -0400 Subject: [PATCH 13/18] Add example to demonstrate noise_learn functionality Some other small modifications were also necessary, e.g., fixing a typo in MCMC.jl, and allowing the user to pass a timestep to the EKIObj. Also, the Cloudy test has been updated to work with the modified code. --- Manifest.toml | 113 ++++-------------------- Project.toml | 3 +- examples/GPEmulator/learn_noise.jl | 134 +++++++++++++++++++++++++++++ examples/GPEmulator/plot_GP.jl | 9 +- src/EKI.jl | 7 +- src/GPEmulator.jl | 6 +- src/MCMC.jl | 8 +- src/Observations.jl | 4 +- test/Cloudy/runtests.jl | 35 ++++---- test/GPEmulator/runtests.jl | 4 +- test/runtests.jl | 2 +- 11 files changed, 195 insertions(+), 130 deletions(-) create mode 100644 examples/GPEmulator/learn_noise.jl diff --git a/Manifest.toml b/Manifest.toml index 6c9b19533..5a54ef3df 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,11 +1,5 @@ # This file is machine-generated - editing it directly is not advised -[[AbstractFFTs]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "051c95d6836228d120f5f4b984dd5aba1624f716" -uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" -version = "0.5.0" - [[AbstractTrees]] deps = ["Markdown"] git-tree-sha1 = "33e450545eaf7699da1a6e755f9ea65f14077a45" @@ -100,9 +94,9 @@ version = "0.8.1" [[ChainRulesCore]] deps = ["MuladdMacro"] -git-tree-sha1 = "4177411bef28d680948562abd25059dceb4d53ed" +git-tree-sha1 = "87e289253a5fc690c4860a41bbd427e16576f716" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "0.9.2" +version = "0.9.3" [[Cloudy]] deps = ["Coverage", "DifferentialEquations", "DocStringExtensions", "ForwardDiff", "HCubature", "LinearAlgebra", "Optim", "PyPlot", "SpecialFunctions", "TaylorSeries", "Test"] @@ -230,9 +224,9 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" [[DiffEqBase]] deps = ["ArrayInterface", "ChainRulesCore", "ConsoleProgressMonitor", "DataStructures", "Distributed", "DocStringExtensions", "FunctionWrappers", "IterativeSolvers", "IteratorInterfaceExtensions", "LabelledArrays", "LinearAlgebra", "Logging", "LoggingExtras", "MuladdMacro", "Parameters", "Printf", "ProgressLogging", "RecipesBase", "RecursiveArrayTools", "RecursiveFactorization", "Requires", "Roots", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "TableTraits", "TerminalLoggers", "TreeViews", "ZygoteRules"] -git-tree-sha1 = "1c977b0c5219e1a4e3ae1ee16b2f29ccc3bd43b7" +git-tree-sha1 = "e7e31a36de80ff2f94bdacd4d956cb6b67fcf0ef" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.40.4" +version = "6.40.7" [[DiffEqCallbacks]] deps = ["DataStructures", "DiffEqBase", "ForwardDiff", "LinearAlgebra", "NLsolve", "OrdinaryDiffEq", "RecipesBase", "RecursiveArrayTools", "StaticArrays"] @@ -310,24 +304,6 @@ git-tree-sha1 = "c5714d9bcdba66389612dc4c47ed827c64112997" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" version = "0.8.2" -[[Documenter]] -deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] -git-tree-sha1 = "395fa1554c69735802bba37d9e7d9586fd44326c" -uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "0.24.11" - -[[ElasticArrays]] -deps = ["Adapt"] -git-tree-sha1 = "b740f7f3fa87a2fa22c93582f2dbb9110c4ece6e" -uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" -version = "1.2.3" - -[[ElasticPDMats]] -deps = ["LinearAlgebra", "MacroTools", "PDMats", "Test"] -git-tree-sha1 = "c6ac4406d9b4a549f7316fc746f7ccd5f1bcd2cd" -uuid = "2904ab23-551e-5aed-883f-487f97af5226" -version = "0.2.1" - [[ExponentialUtilities]] deps = ["LinearAlgebra", "Printf", "Requires", "SparseArrays"] git-tree-sha1 = "91f7498b66205431fe3e35833cda97a22b1ab6a5" @@ -346,29 +322,11 @@ git-tree-sha1 = "0fa07f43e5609ea54848b82b4bb330b250e9645b" uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" version = "4.1.0+3" -[[FFTW]] -deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"] -git-tree-sha1 = "14536c95939aadcee44014728a459d2fe3ca9acf" -uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -version = "1.2.2" - -[[FFTW_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "6c975cd606128d45d1df432fb812d6eb10fee00b" -uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" -version = "3.3.9+5" - -[[FastGaussQuadrature]] -deps = ["LinearAlgebra", "SpecialFunctions"] -git-tree-sha1 = "c139e3f4c75dc489a493627c7ee44befc177420f" -uuid = "442a2c76-b920-505d-bb47-c5924d526838" -version = "0.4.2" - [[FillArrays]] deps = ["LinearAlgebra", "Random", "SparseArrays"] -git-tree-sha1 = "4783bbbeade37f2a8bd82af6c112510fde78e532" +git-tree-sha1 = "be4180bdb27a11188d694ee3773122f4921f1a62" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "0.8.12" +version = "0.8.13" [[FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] @@ -377,9 +335,10 @@ uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" version = "2.4.1" [[FixedPointNumbers]] -git-tree-sha1 = "8fb797c37a3b7ced4327a05ac4ca0dd6a4f1ba92" +deps = ["Statistics"] +git-tree-sha1 = "266baee2e9d875cb7a3bfdcc6cab553c543ff8ab" uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" -version = "0.8.1" +version = "0.8.2" [[Formatting]] deps = ["Printf"] @@ -420,12 +379,6 @@ git-tree-sha1 = "247adbd2b33c0c4b42efa20d1e807acf6312145f" uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" version = "0.50.1" -[[GaussianProcesses]] -deps = ["Distances", "Distributions", "Documenter", "ElasticArrays", "ElasticPDMats", "FastGaussQuadrature", "ForwardDiff", "LinearAlgebra", "Optim", "PDMats", "Printf", "ProgressMeter", "Random", "RecipesBase", "ScikitLearnBase", "SpecialFunctions", "StaticArrays", "Statistics", "StatsFuns", "Zygote"] -git-tree-sha1 = "9b33cd8be0fd2bc59c46911def57331a7d5f6e84" -uuid = "891a1506-143c-57d2-908e-e1f8e92e6de9" -version = "0.12.1" - [[GeneralizedGenerated]] deps = ["CanonicalTraits", "DataStructures", "JuliaVariables", "MLStyle"] git-tree-sha1 = "9a917449bebc5a241d23e13d36ea62c77129ce6e" @@ -462,12 +415,6 @@ git-tree-sha1 = "eca61b35cdd8cd2fcc5eec1eda766424a995b02f" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" version = "0.8.16" -[[IRTools]] -deps = ["InteractiveUtils", "MacroTools", "Test"] -git-tree-sha1 = "90ee39f9beaaa186e4968417ea2b8ed5673c91c0" -uuid = "7869d1d1-7146-5819-86e3-90919afe41df" -version = "0.3.3" - [[Inflate]] git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c" uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" @@ -479,12 +426,6 @@ git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8" uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f" version = "0.5.0" -[[IntelOpenMP_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "fb8e1c7a5594ba56f9011310790e03b5384998d6" -uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" -version = "2018.0.3+0" - [[InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" @@ -596,9 +537,9 @@ version = "0.4.1" [[LoopVectorization]] deps = ["DocStringExtensions", "LinearAlgebra", "OffsetArrays", "SIMDPirates", "SLEEFPirates", "UnPack", "VectorizationBase"] -git-tree-sha1 = "6f639e9cb9273c500f18779e4805646e6ac57d21" +git-tree-sha1 = "b595e15d20e45d2eb36c6b4462d2a34143872a45" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.8.14" +version = "0.8.15" [[METIS_jll]] deps = ["Libdl", "Pkg"] @@ -606,12 +547,6 @@ git-tree-sha1 = "3f52ed323683398498ef163a45ce998f1ceca363" uuid = "d00139f3-1899-568f-a2f0-47f597d42d70" version = "5.1.0+4" -[[MKL_jll]] -deps = ["IntelOpenMP_jll", "Libdl", "Pkg"] -git-tree-sha1 = "0ce9a7fa68c70cf83c49d05d2c04d91b47404b08" -uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -version = "2020.1.216+0" - [[MLStyle]] git-tree-sha1 = "67f9a88611bc79f992aa705d9bbc833a2547dec7" uuid = "d8e11817-5142-5d16-987a-aa16d5891078" @@ -682,16 +617,10 @@ git-tree-sha1 = "ea172c86745810136d744fc67650d2e2de669c4f" uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" version = "4.4.0" -[[NNlib]] -deps = ["BinaryProvider", "Libdl", "LinearAlgebra", "Requires", "Statistics"] -git-tree-sha1 = "d9f196d911f55aeaff11b11f681b135980783824" -uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" -version = "0.6.6" - [[NaNMath]] -git-tree-sha1 = "928b8ca9b2791081dc71a51c55347c27c618760f" +git-tree-sha1 = "c84c576296d0e2fbb3fc134d3e09086b3ea617cd" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "0.3.3" +version = "0.3.4" [[NameResolution]] deps = ["DataStructures", "PrettyPrint"] @@ -752,10 +681,10 @@ uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" version = "0.9.12" [[ParameterizedFunctions]] -deps = ["DataStructures", "DiffEqBase", "Latexify", "LinearAlgebra", "ModelingToolkit"] -git-tree-sha1 = "86ca1fa2f37f20945a742d50ed12e314b11dfa62" +deps = ["DataStructures", "DiffEqBase", "Latexify", "LinearAlgebra", "ModelingToolkit", "Reexport"] +git-tree-sha1 = "458f179b2bcb8071a2dbcf92cafd2c1cbd44e8e1" uuid = "65888b18-ceab-5e60-b2b9-181511a3b968" -version = "5.3.0" +version = "5.4.0" [[Parameters]] deps = ["OrderedCollections", "UnPack"] @@ -826,9 +755,9 @@ version = "0.1.3" [[ProgressMeter]] deps = ["Distributed", "Printf"] -git-tree-sha1 = "3e1784c27847bba115815d4d4e668b99873985e5" +git-tree-sha1 = "2de4cddc0ceeddafb6b143b5b6cd9c659b64507c" uuid = "92933f4c-e287-5a05-a399-4b506db050ca" -version = "1.3.1" +version = "1.3.2" [[PyCall]] deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"] @@ -1141,12 +1070,6 @@ git-tree-sha1 = "622d8b6dc0c7e8029f17127703de9819134d1b71" uuid = "83775a58-1f1d-513f-b197-d71354ab007a" version = "1.2.11+14" -[[Zygote]] -deps = ["DiffRules", "FFTW", "FillArrays", "ForwardDiff", "IRTools", "InteractiveUtils", "LinearAlgebra", "MacroTools", "NNlib", "NaNMath", "Random", "Requires", "SpecialFunctions", "Statistics", "ZygoteRules"] -git-tree-sha1 = "f8329b595c465caf3ca87c4f744e6041a4983e43" -uuid = "e88e6eb3-aa80-5325-afca-941959d7151f" -version = "0.4.8" - [[ZygoteRules]] deps = ["MacroTools"] git-tree-sha1 = "b3b4882cc9accf6731a08cc39543fbc6b669dca8" diff --git a/Project.toml b/Project.toml index 2a5f22514..4d3a771d0 100644 --- a/Project.toml +++ b/Project.toml @@ -9,14 +9,13 @@ Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Optim = "429524aa-4258-5aef-a3af-852621145aeb" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/examples/GPEmulator/learn_noise.jl b/examples/GPEmulator/learn_noise.jl new file mode 100644 index 000000000..39d2b82c2 --- /dev/null +++ b/examples/GPEmulator/learn_noise.jl @@ -0,0 +1,134 @@ +# Import modules +using Random +using GaussianProcesses +using Distributions +using Statistics +using LinearAlgebra +using CalibrateEmulateSample.GPEmulator + +############################################################################### +# # +# This examples shows how to fit a Gaussian Process (GP) regression model # +# using GPEmulator and demonstrates how the GPObj's built-in Singular # +# Value Decomposition (SVD) decorrelates the training data and maps into # +# a space where the covariance of the observational noise is the identity. # +# # +# When training a GP, the hyperparameters of its kernel (or covariance # +# function) are optimized with respect to the given input-output pairs. # +# When training a separate GP for each output component (as is usually # +# done for multidimensional output, rather than fitting one multi-output # +# model), one assumes that their covariance functions are independent. # +# The SVD transformation ensures that this is a valid assumption. # +# # +# The decorrelation by the SVD is done automatically when instantiating # +# a `GPObj`, but the user can choose to have the `predict` function # +# return its predictions in the original space (by setting # +# transform_to_real=true) # +# # +# Training points: {(x_i, y_i)} (i=1,...,n), where x_i ∈ ℝ ² and y_i ∈ ℝ ² # +# The y_i are assumed to be related to the x_i by: # +# y_i = G(x_i) + η, where G is a map from ℝ ² to ℝ ², and η is noise drawn # +# from a 2D normal distribution with known mean μ and covariance Σ # +# Two Gaussian Process models are fit, one that predicts y_i[1] from x_i # +# and one that predicts y_i[2] from x_i. # +# # +# The GPEmulator module can be used as a standalone module to fit # +# Gaussian Process regression models, but it was originally designed as # +# the "Emulate" step of the "Calibrate-Emulate-Sample" framework # +# developed at CliMA. # +# Further reading on Calibrate-Emulate-Sample: # +# Cleary et al., 2020 # +# https://arxiv.org/abs/2001.03689 # +# # +############################################################################### + +rng_seed = 41 +Random.seed!(rng_seed) + +gppackage = GPJL() +pred_type = YType() + +# Generate training data (x-y pairs, where x ∈ ℝ ᵖ, y ∈ ℝ ᵈ) +# x = [x1, x2]: inputs/predictors/features/parameters +# y = [y1, y2]: outputs/predictands/targets/observables +# The observables y are related to the parameters x by: +# y = G(x1, x2) + η, +# where G(x1, x2) := [sin(x1) + cos(x2), sin(x1) - cos(x2)], and η ~ N(0, Σ) +n = 50 # number of training points +p = 2 # input dim +d = 2 # output dim + +X = 2.0 * π * rand(n, p) + +# G(x1, x2) +g1x = sin.(X[:, 1]) .+ cos.(X[:, 2]) +g2x = sin.(X[:, 1]) .- cos.(X[:, 2]) +gx = [g1x g2x] + +# Add noise η +μ = zeros(d) +Σ = 0.1 * [[0.8, 0.2] [0.2, 0.5]] # d x d +noise_samples = rand(MvNormal(μ, Σ), n)' + +# y = G(x) + η +Y = gx .+ noise_samples +input_cov = cov(noise_samples, dims=1) + +# Fit 2D Gaussian Process regression model +# (To be precise: We fit two models, one that predicts y1 from x1 and x2, and +# one that predicts y2 from x1 and x2) +# Setting GPkernel=nothing leads to the creation of a default squared +# exponential kernel. +# Setting noise_learn=true leads to the addition of white noise to the kernel, +# whose hyperparameter -- the signal variance, or "noise" -- will be learned +# in the training phase via an optimization procedure. +# Because of the SVD transformation applied to the output, we expect the +# learned noise to be close to 1. +gpobj1 = GPObj(X, Y, input_cov, gppackage, GPkernel=nothing, + normalized=true, noise_learn=true, prediction_type=pred_type) + +println("\n-----------") +println("Results of training Gaussian Process models with noise_learn=true") +println("-----------") + +println("\nKernel of the GP trained to predict y1 from x=(x1, x2):") +println(gpobj1.models[1].kernel) +println("Learned noise parameter, σ_1:") +learned_σ_1 = sqrt(gpobj1.models[1].kernel.kright.σ2) +println("σ_1 = $learned_σ_1") +# Check if the learned noise is approximately 1 +@assert(isapprox(learned_σ_1, 1.0; atol=0.07)) + +println("\nKernel of the GP trained to predict y2 from x=(x1, x2):") +println(gpobj1.models[2].kernel) +println("Learned noise parameter, σ_2:") +learned_σ_2 = sqrt(gpobj1.models[2].kernel.kright.σ2) +println("σ_2 = $learned_σ_2") +# Check if the learned noise is approximately 1 +@assert(isapprox(learned_σ_2, 1.0; atol=0.07)) +println("------------------------------------------------------------------\n") + +# For comparison: When noise_learn is set to false, the observational noise +# is set to 1.0 and is not learned/optimized during the training. But thanks +# to the SVD, 1.0 is the correct value to use. +gpobj2 = GPObj(X, Y, input_cov, gppackage, GPkernel=nothing, + normalized=true, noise_learn=false, prediction_type=pred_type) + +println("\n-----------") +println("Results of training Gaussian Process models with noise_learn=false") +println("-----------") + +println("\nKernel of the GP trained to predict y1 from x=(x1, x2):") +# Note: In contrast to the kernels of the gpobj1 models, these ones do not +# have a white noise ("Noise") kernel component +println(gpobj2.models[1].kernel) +# logNoise is given as log(sqrt(noise)) +obs_noise_1 = exp(gpobj2.models[1].logNoise.value^2) +println("Observational noise: $obs_noise_1") + +println("\nKernel of the GP trained to predict y1 from x=(x1, x2):") +println(gpobj2.models[2].kernel) +# logNoise is given as log(sqrt(noise)) +obs_noise_2 = exp(gpobj2.models[2].logNoise.value^2) +println("Observational noise: $obs_noise_2") +println("------------------------------------------------------------------") diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl index e1344623a..0752f4e7a 100644 --- a/examples/GPEmulator/plot_GP.jl +++ b/examples/GPEmulator/plot_GP.jl @@ -22,8 +22,9 @@ using CalibrateEmulateSample.GPEmulator # and one that predicts y_i[2] from x_i. Fitting separate models for # # y_i[1] and y_i[2] requires the outputs to be independent - since this # # cannot be assumed to be the case a priori, the training data are # -# decorrelated by perfoming an SVD decomposition on the noise # -# covariance Σ, and each model is then trained in the decorrelated space. # +# decorrelated by perfoming an Singulat Value Decomposition (SVD) on the # +# noise covariance Σ, and each model is then trained in the decorrelated # +# space. # # # # The decorrelation is done automatically when instantiating a `GPObj`, # # but the user can choose to have the `predict` function return its # @@ -77,7 +78,7 @@ noise_samples = rand(MvNormal(μ, Σ), n)' # y = G(x) + η Y = gx .+ noise_samples -input_cov = cov(noise_samples, dims=1) +data_cov = cov(noise_samples, dims=1) # Fit 2D Gaussian Process regression model # (To be precise: We fit two models, one that predicts y1 from x1 and x2, @@ -86,7 +87,7 @@ input_cov = cov(noise_samples, dims=1) # exponential kernel. # Setting noise_learn=true leads to the addition of white noise to the # kernel -gpobj = GPObj(X, Y, input_cov, gppackage, GPkernel=nothing, +gpobj = GPObj(X, Y, data_cov, gppackage, GPkernel=nothing, normalized=true, noise_learn=true, prediction_type=pred_type) # Plot mean and variance of the predicted observables y1 and y2 diff --git a/src/EKI.jl b/src/EKI.jl index b440dc408..455f10ceb 100644 --- a/src/EKI.jl +++ b/src/EKI.jl @@ -44,7 +44,8 @@ end function EKIObj(parameters::Array{FT, 2}, parameter_names::Vector{String}, t_mean, - t_cov::Array{FT, 2}) where {FT<:AbstractFloat} + t_cov::Array{FT, 2}; + Δt=FT(1)) where {FT<:AbstractFloat} # ensemble size N_ens = size(parameters)[1] @@ -57,7 +58,7 @@ function EKIObj(parameters::Array{FT, 2}, # error store err = FT[] # timestep store - Δt = FT[] + Δt = Array([Δt]) EKIObj{FT,IT}(u, parameter_names, t_mean, t_cov, N_ens, g, err, Δt) end @@ -95,7 +96,7 @@ function compute_error(eki) end -function update_ensemble!(eki::EKIObj{FT}, g; cov_threshold::FT=0.01, Δt_new = nothing) where {FT} +function update_ensemble!(eki::EKIObj{FT}, g; cov_threshold::FT=0.01, Δt_new=nothing) where {FT} # u: N_ens x N_params u = eki.u[end] cov_init = cov(eki.u[end], dims=1) diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 37ffd7983..7ecf610aa 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -124,7 +124,7 @@ function GPObj( GPinputs = convert(Array, inputs') end - # Transform the outputs (transformed_data will also come out transposed from + # Transform the outputs (transformed_data will come out transposed from # this transformation, i.e., with size output_dim x N_samples. This is also # the size that will be required when passing the data to GPE(), so we leave # transformed_data in that size. @@ -152,7 +152,7 @@ function GPObj( end if noise_learn - # Add white noise to kernel + # Add white noise to kernel white_logstd = log(1.0) white = Noise(white_logstd) kern = kern + white @@ -172,7 +172,9 @@ function GPObj( for i in 1:N_models # Make a copy of the kernel (because it gets altered in every # iteration) + println("kernel in GPObj") GPkernel_i = deepcopy(kern) + println(GPkernel_i) GPdata_i = transformed_data[i, :] # GPE() arguments: # GPinputs: (input_dim x N_samples) diff --git a/src/MCMC.jl b/src/MCMC.jl index e8913772c..1d827f311 100644 --- a/src/MCMC.jl +++ b/src/MCMC.jl @@ -82,12 +82,12 @@ function MCMCObj( param_init_copy = deepcopy(param_init) - #We need to transform obs_sample into the correct space ( + # We need to transform obs_sample into the correct space if svdflag println("Applying SVD to decorrelating outputs, if not required set svdflag=false") - decomposition=svd(truth_cov)#svd.U * svd.S * svd.Vt (can also get V) - sqrt_singular_values_inv=Diagonal(1.0 ./ sqrt.(decomposition.S)) #diagonal matrix of 1/eigenvalues - obs_sample=sqrt_singular_values_inv*decomposition.Vt*obs_sample + decomposition = svd(obs_cov) # svd.U * svd.S * svd.Vt (can also get V) + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) # diagonal matrix of 1/eigenvalues + obs_sample = sqrt_singular_values_inv * decomposition.Vt * obs_sample else println("Assuming independent outputs.") end diff --git a/src/Observations.jl b/src/Observations.jl index cc8b63675..c1ba2b2ed 100644 --- a/src/Observations.jl +++ b/src/Observations.jl @@ -42,7 +42,7 @@ function Obs(samples::Vector{Vector{FT}}, else temp = convert(Array, reshape(hcat(samples...)', N_samples, :)) samplemean = vec(mean(temp, dims=1)) - cov = Statistics.cov(temp) + cov = Statistics.cov(temp .- samplemean) end Obs(samples, cov, samplemean, data_names) end @@ -57,10 +57,10 @@ function Obs(samples::Array{FT}, samplemean = vec(samples) samples = vec([vec(samples)]) else - cov = Statistics.cov(samples) # convert matrix of samples to a vector of vectors N_samples = size(samples, 1) samplemean = vec(mean(samples, dims=1)) + cov = Statistics.cov(samples .- samplemean) samples = [samples[i, :] for i in 1:N_samples] end Obs(samples, cov, samplemean, data_names) diff --git a/test/Cloudy/runtests.jl b/test/Cloudy/runtests.jl index 7c947f395..3d218da2e 100644 --- a/test/Cloudy/runtests.jl +++ b/test/Cloudy/runtests.jl @@ -19,6 +19,7 @@ using CalibrateEmulateSample.MCMC using CalibrateEmulateSample.Observations using CalibrateEmulateSample.GModel using CalibrateEmulateSample.Utilities +using CalibrateEmulateSample.Priors ### @@ -84,23 +85,25 @@ tspan = (0., 0.5) ### ### Generate (artificial) truth samples +### Note: The observables y are related to the parameters θ by: +### y = G(x1, x2) + η ### g_settings_true = GModel.GSettings(kernel, dist_true, moments, tspan) -yt = GModel.run_G(params_true, g_settings_true, PDistributions.update_params, +gt = GModel.run_G(params_true, g_settings_true, PDistributions.update_params, PDistributions.moment, Cloudy.Sources.get_int_coalescence) n_samples = 100 -samples = zeros(n_samples, length(yt)) +yt = zeros(n_samples, length(gt)) noise_level = 0.05 -Γy = noise_level^2 * convert(Array, Diagonal(yt)) -μ = zeros(length(yt)) +Γy = noise_level * convert(Array, Diagonal(gt)) +μ = zeros(length(gt)) # Add noise for i in 1:n_samples - samples[i, :] = yt + noise_level^2 * rand(MvNormal(μ, Γy)) + yt[i, :] = gt .+ rand(MvNormal(μ, Γy)) end -truth = Observations.Obs(samples, Γy, data_names) +truth = Observations.Obs(yt, Γy, data_names) ### @@ -114,7 +117,7 @@ N_ens = 50 # number of ensemble members N_iter = 5 # number of EKI iterations # initial parameters: N_ens x N_params initial_params = EKI.construct_initial_ensemble(N_ens, priors; rng_seed=6) -ekiobj = EKI.EKIObj(initial_params, param_names, truth.mean, truth.cov) +ekiobj = EKI.EKIObj(initial_params, param_names, truth.mean, truth.cov, Δt=1.0) # Initialize a ParticleDistribution with dummy parameters. The parameters # will then be set in run_G_ensemble @@ -159,17 +162,19 @@ len2 = zeros(3) kern2 = Mat52Ard(len2, 0.0) # regularize with white noise white = Noise(log(2.0)) -# construct kernel +# # construct kernel GPkernel = kern1 + kern2 + white # Get training points u_tp, g_tp = Utilities.extract_GP_tp(ekiobj, N_iter) normalized = true -gpobj = GPEmulator.GPObj(u_tp, g_tp, gppackage; GPkernel=GPkernel, - normalized=normalized, prediction_type=pred_type) +gpobj = GPEmulator.GPObj(u_tp, g_tp, Γy, gppackage; GPkernel=GPkernel, + normalized=normalized, noise_learn=false, + prediction_type=pred_type) # Check how well the Gaussian Process regression predicts on the # true parameters -y_mean, y_var = GPEmulator.predict(gpobj, reshape(log.(params_true), 1, :)) +y_mean, y_var = GPEmulator.predict(gpobj, reshape(log.(params_true), 1, :), + transform_to_real=true) println("GP prediction on true parameters: ") println(vec(y_mean)) @@ -193,8 +198,8 @@ burnin = 0 step = 0.1 # first guess max_iter = 5000 yt_sample = truth.mean -mcmc_test = MCMC.MCMCObj(yt_sample, truth.cov, priors, step, u0, max_iter, - mcmc_alg, burnin) +mcmc_test = MCMC.MCMCObj(yt_sample, Γy, priors, step, u0, max_iter, + mcmc_alg, burnin, svdflag=true) new_step = MCMC.find_mcmc_step!(mcmc_test, gpobj) # Now begin the actual MCMC @@ -205,8 +210,8 @@ u0 = vec(mean(u_tp, dims=1)) burnin = 1000 max_iter = 100000 -mcmc = MCMC.MCMCObj(yt_sample, truth.cov, priors, - new_step, u0, max_iter, mcmc_alg, burnin) +mcmc = MCMC.MCMCObj(yt_sample, Γy, priors, new_step, u0, max_iter, + mcmc_alg, burnin, svdflag=true) MCMC.sample_posterior!(mcmc, gpobj, max_iter) posterior = MCMC.get_posterior(mcmc) diff --git a/test/GPEmulator/runtests.jl b/test/GPEmulator/runtests.jl index 4abb7f238..143f76466 100644 --- a/test/GPEmulator/runtests.jl +++ b/test/GPEmulator/runtests.jl @@ -28,8 +28,8 @@ using CalibrateEmulateSample.GPEmulator GPkernel = kern + white # Fit Gaussian Process Regression model - gpobj = GPObj(x, y, gppackage; GPkernel=GPkernel, normalized=false, - prediction_type=pred_type) + gpobj = GPObj(x, y, input_cov, gppackage; GPkernel=GPkernel, normalized=true, + noise_learn=true, prediction_type=pred_type) new_inputs = [0.0, π/2, π, 3*π/2, 2*π] μ, σ² = GPEmulator.predict(gpobj, new_inputs) diff --git a/test/runtests.jl b/test/runtests.jl index 717af6e7a..8571cbc82 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,7 +14,7 @@ for submodule in [ #"GPR", #"EKS", #"Histograms", - #"Cloudy", + "Cloudy", ] println("Starting tests for $submodule") From 5c7e89c5afaa4ce91dfd21efdf6a45771068180a Mon Sep 17 00:00:00 2001 From: Melanie Bieli Date: Thu, 16 Jul 2020 00:46:00 -0400 Subject: [PATCH 14/18] Get all tests to work and major clean-up --- Manifest.toml | 187 ++++++++++++++-- Project.toml | 3 + examples/GPEmulator/learn_noise.jl | 3 +- examples/GPEmulator/plot_GP.jl | 3 +- src/EKI.jl | 4 +- src/GPEmulator.jl | 343 +++++++++++++++++------------ src/MCMC.jl | 43 ++-- src/Observations.jl | 135 ++++++++---- test/Cloudy/runtests.jl | 12 +- test/EKI/runtests.jl | 2 +- test/GPEmulator/runtests.jl | 148 ++++++++++--- test/MCMC/runtests.jl | 40 ++-- test/Observations/runtests.jl | 44 +++- 13 files changed, 677 insertions(+), 290 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 5a54ef3df..53f6ea3e1 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,5 +1,11 @@ # This file is machine-generated - editing it directly is not advised +[[AbstractFFTs]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "051c95d6836228d120f5f4b984dd5aba1624f716" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "0.5.0" + [[AbstractTrees]] deps = ["Markdown"] git-tree-sha1 = "33e450545eaf7699da1a6e755f9ea65f14077a45" @@ -38,9 +44,15 @@ version = "2.9.1" [[ArrayLayouts]] deps = ["FillArrays", "LinearAlgebra"] -git-tree-sha1 = "a3254b3780a3544838ca0b7e23b1e9b06eb71bd8" +git-tree-sha1 = "6f6e33efac70fc24c1f2a654a090b7af01690ffe" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" -version = "0.3.5" +version = "0.3.7" + +[[AxisAlgorithms]] +deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] +git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9" +uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950" +version = "1.0.0" [[BandedMatrices]] deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"] @@ -104,6 +116,12 @@ git-tree-sha1 = "7cb3ce76f5a63d17eb358b6030f7664931c75857" uuid = "9e3b23bb-e7cc-4b94-886c-65de2234ba87" version = "0.1.0" +[[Clustering]] +deps = ["Distances", "LinearAlgebra", "NearestNeighbors", "Printf", "SparseArrays", "Statistics", "StatsBase"] +git-tree-sha1 = "b11c8d607af357776a046889a7c32567d05f1319" +uuid = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" +version = "0.14.1" + [[ColorSchemes]] deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"] git-tree-sha1 = "7a15e3690529fd1042f0ab954dff7445b1efc8a5" @@ -208,6 +226,12 @@ git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464" version = "1.0.0" +[[DataValues]] +deps = ["DataValueInterfaces", "Dates"] +git-tree-sha1 = "d88a19299eba280a6d062e135a43f00323ae70bf" +uuid = "e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5" +version = "0.4.13" + [[Dates]] deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" @@ -230,9 +254,9 @@ version = "6.40.7" [[DiffEqCallbacks]] deps = ["DataStructures", "DiffEqBase", "ForwardDiff", "LinearAlgebra", "NLsolve", "OrdinaryDiffEq", "RecipesBase", "RecursiveArrayTools", "StaticArrays"] -git-tree-sha1 = "c0619bba19db9b033d011b91e30663e4baf74184" +git-tree-sha1 = "1b7d2b06490bd917c3c5e7bf65e6490537bb1652" uuid = "459566f4-90b8-5000-8ac3-15dfb0a30def" -version = "2.13.3" +version = "2.13.4" [[DiffEqFinancial]] deps = ["DiffEqBase", "DiffEqNoiseProcess", "LinearAlgebra", "Markdown", "RandomNumbers"] @@ -304,6 +328,24 @@ git-tree-sha1 = "c5714d9bcdba66389612dc4c47ed827c64112997" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" version = "0.8.2" +[[Documenter]] +deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] +git-tree-sha1 = "395fa1554c69735802bba37d9e7d9586fd44326c" +uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +version = "0.24.11" + +[[ElasticArrays]] +deps = ["Adapt"] +git-tree-sha1 = "b740f7f3fa87a2fa22c93582f2dbb9110c4ece6e" +uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" +version = "1.2.3" + +[[ElasticPDMats]] +deps = ["LinearAlgebra", "MacroTools", "PDMats", "Test"] +git-tree-sha1 = "c6ac4406d9b4a549f7316fc746f7ccd5f1bcd2cd" +uuid = "2904ab23-551e-5aed-883f-487f97af5226" +version = "0.2.1" + [[ExponentialUtilities]] deps = ["LinearAlgebra", "Printf", "Requires", "SparseArrays"] git-tree-sha1 = "91f7498b66205431fe3e35833cda97a22b1ab6a5" @@ -322,6 +364,24 @@ git-tree-sha1 = "0fa07f43e5609ea54848b82b4bb330b250e9645b" uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" version = "4.1.0+3" +[[FFTW]] +deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"] +git-tree-sha1 = "14536c95939aadcee44014728a459d2fe3ca9acf" +uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +version = "1.2.2" + +[[FFTW_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "6c975cd606128d45d1df432fb812d6eb10fee00b" +uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" +version = "3.3.9+5" + +[[FastGaussQuadrature]] +deps = ["LinearAlgebra", "SpecialFunctions"] +git-tree-sha1 = "c139e3f4c75dc489a493627c7ee44befc177420f" +uuid = "442a2c76-b920-505d-bb47-c5924d526838" +version = "0.4.2" + [[FillArrays]] deps = ["LinearAlgebra", "Random", "SparseArrays"] git-tree-sha1 = "be4180bdb27a11188d694ee3773122f4921f1a62" @@ -379,6 +439,12 @@ git-tree-sha1 = "247adbd2b33c0c4b42efa20d1e807acf6312145f" uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" version = "0.50.1" +[[GaussianProcesses]] +deps = ["Distances", "Distributions", "Documenter", "ElasticArrays", "ElasticPDMats", "FastGaussQuadrature", "ForwardDiff", "LinearAlgebra", "Optim", "PDMats", "Printf", "ProgressMeter", "Random", "RecipesBase", "ScikitLearnBase", "SpecialFunctions", "StaticArrays", "Statistics", "StatsFuns", "Zygote"] +git-tree-sha1 = "9b33cd8be0fd2bc59c46911def57331a7d5f6e84" +uuid = "891a1506-143c-57d2-908e-e1f8e92e6de9" +version = "0.12.1" + [[GeneralizedGenerated]] deps = ["CanonicalTraits", "DataStructures", "JuliaVariables", "MLStyle"] git-tree-sha1 = "9a917449bebc5a241d23e13d36ea62c77129ce6e" @@ -415,6 +481,12 @@ git-tree-sha1 = "eca61b35cdd8cd2fcc5eec1eda766424a995b02f" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" version = "0.8.16" +[[IRTools]] +deps = ["InteractiveUtils", "MacroTools", "Test"] +git-tree-sha1 = "90ee39f9beaaa186e4968417ea2b8ed5673c91c0" +uuid = "7869d1d1-7146-5819-86e3-90919afe41df" +version = "0.3.3" + [[Inflate]] git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c" uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" @@ -426,10 +498,22 @@ git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8" uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f" version = "0.5.0" +[[IntelOpenMP_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "fb8e1c7a5594ba56f9011310790e03b5384998d6" +uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" +version = "2018.0.3+0" + [[InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +[[Interpolations]] +deps = ["AxisAlgorithms", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] +git-tree-sha1 = "2b7d4e9be8b74f03115e64cf36ed2f48ae83d946" +uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +version = "0.12.10" + [[InvertedIndices]] deps = ["Test"] git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc" @@ -464,18 +548,18 @@ git-tree-sha1 = "8868479ff35ab98588ed0a529a9c2a4f8bda3bd6" uuid = "b14d175d-62b4-44ba-8fb7-3064adc8c3ec" version = "0.2.0" +[[KernelDensity]] +deps = ["Distributions", "FFTW", "Interpolations", "Optim", "StatsBase", "Test"] +git-tree-sha1 = "c1048817fe5711f699abc8fabd47b1ac6ba4db04" +uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" +version = "0.5.1" + [[LAME_jll]] deps = ["Libdl", "Pkg"] git-tree-sha1 = "221cc8998b9060677448cbb6375f00032554c4fd" uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d" version = "3.100.0+1" -[[LLVM]] -deps = ["CEnum", "Libdl", "Printf", "Unicode"] -git-tree-sha1 = "a662366a5d485dee882077e8da3e1a95a86d097f" -uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "2.0.0" - [[LaTeXStrings]] git-tree-sha1 = "de44b395389b84fd681394d4e8d39ef14e3a2ea8" uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" @@ -531,15 +615,16 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[LoggingExtras]] -git-tree-sha1 = "b60616c70eff0cc2c0831b6aace75940aeb0939d" +deps = ["Dates"] +git-tree-sha1 = "03289aba73c0abc25ff0229bed60f2a4129cd15c" uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" -version = "0.4.1" +version = "0.4.2" [[LoopVectorization]] deps = ["DocStringExtensions", "LinearAlgebra", "OffsetArrays", "SIMDPirates", "SLEEFPirates", "UnPack", "VectorizationBase"] -git-tree-sha1 = "b595e15d20e45d2eb36c6b4462d2a34143872a45" +git-tree-sha1 = "4c002de66221639174e081ec74c156e9013e4afa" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.8.15" +version = "0.8.17" [[METIS_jll]] deps = ["Libdl", "Pkg"] @@ -547,6 +632,12 @@ git-tree-sha1 = "3f52ed323683398498ef163a45ce998f1ceca363" uuid = "d00139f3-1899-568f-a2f0-47f597d42d70" version = "5.1.0+4" +[[MKL_jll]] +deps = ["IntelOpenMP_jll", "Libdl", "Pkg"] +git-tree-sha1 = "0ce9a7fa68c70cf83c49d05d2c04d91b47404b08" +uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +version = "2020.1.216+0" + [[MLStyle]] git-tree-sha1 = "67f9a88611bc79f992aa705d9bbc833a2547dec7" uuid = "d8e11817-5142-5d16-987a-aa16d5891078" @@ -605,6 +696,12 @@ git-tree-sha1 = "258f3be6770fe77be8870727ba9803e236c685b8" uuid = "f9640e96-87f6-5992-9c3b-0743c6a49ffa" version = "1.8.1" +[[MultivariateStats]] +deps = ["Arpack", "LinearAlgebra", "SparseArrays", "Statistics", "StatsBase"] +git-tree-sha1 = "352fae519b447bf52e6de627b89f448bcd469e4e" +uuid = "6f286f6a-111f-5878-ab1e-185364afe411" +version = "0.7.0" + [[NLSolversBase]] deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] git-tree-sha1 = "7c4e66c47848562003250f28b579c584e55becc0" @@ -617,6 +714,12 @@ git-tree-sha1 = "ea172c86745810136d744fc67650d2e2de669c4f" uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" version = "4.4.0" +[[NNlib]] +deps = ["BinaryProvider", "Libdl", "LinearAlgebra", "Requires", "Statistics"] +git-tree-sha1 = "d9f196d911f55aeaff11b11f681b135980783824" +uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" +version = "0.6.6" + [[NaNMath]] git-tree-sha1 = "c84c576296d0e2fbb3fc134d3e09086b3ea617cd" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" @@ -628,6 +731,17 @@ git-tree-sha1 = "f4119274d5a410c64a1d9f546312bb6ae54d41c0" uuid = "71a1bf82-56d0-4bbc-8a3c-48b961074391" version = "0.1.3" +[[NearestNeighbors]] +deps = ["Distances", "StaticArrays"] +git-tree-sha1 = "93107e3cdada73d63245ed8170dcae680f0c8fd8" +uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce" +version = "0.4.6" + +[[Observables]] +git-tree-sha1 = "11832878355305984235a2e90d0e3737383c634c" +uuid = "510215fc-4207-5dde-b226-833fc4488ee2" +version = "0.3.1" + [[OffsetArrays]] git-tree-sha1 = "4ba4cd84c88df8340da1c3e2d8dcb9d18dd1b53b" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" @@ -791,6 +905,11 @@ git-tree-sha1 = "441e6fc35597524ada7f85e13df1f4e10137d16f" uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143" version = "1.4.0" +[[Ratios]] +git-tree-sha1 = "37d210f612d70f3f7d57d488cb3b6eff56ad4e41" +uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439" +version = "0.4.0" + [[RecipesBase]] git-tree-sha1 = "54f8ceb165a0f6d083f0d12cb4996f5367c6edbc" uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" @@ -855,9 +974,9 @@ uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" [[SIMDPirates]] deps = ["VectorizationBase"] -git-tree-sha1 = "dae629e96c1819d77882256e6cb29736f493bc30" +git-tree-sha1 = "cbe8797ac354d0b1dfe75d983429938db834b706" uuid = "21efa798-c60a-11e8-04d3-e1a92915a26a" -version = "0.8.13" +version = "0.8.16" [[SLEEFPirates]] deps = ["Libdl", "SIMDPirates", "VectorizationBase"] @@ -949,6 +1068,12 @@ git-tree-sha1 = "04a5a8e6ab87966b43f247920eab053fd5fdc925" uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" version = "0.9.5" +[[StatsPlots]] +deps = ["Clustering", "DataStructures", "DataValues", "Distributions", "Interpolations", "KernelDensity", "MultivariateStats", "Observables", "Plots", "RecipesBase", "RecipesPipeline", "Reexport", "StatsBase", "TableOperations", "Tables", "Widgets"] +git-tree-sha1 = "b9b7fff81f573465fcac4685df1497d968537a9e" +uuid = "f3b207a7-027a-5e70-b257-86293d7955fd" +version = "0.14.6" + [[SteadyStateDiffEq]] deps = ["DiffEqBase", "DiffEqCallbacks", "LinearAlgebra", "NLsolve", "Reexport"] git-tree-sha1 = "75f258513b7ef8b235876f4cf146577ffd545094" @@ -989,6 +1114,12 @@ git-tree-sha1 = "4a4ae2ebefa34779552dbe87683c70c856624ec8" uuid = "fb77eaff-e24c-56d4-86b1-d163f2edb164" version = "5.2.0+0" +[[TableOperations]] +deps = ["Tables", "Test"] +git-tree-sha1 = "208630a14884abd110a8f8008b0882f0d0f5632c" +uuid = "ab02a1b2-a7df-11e8-156e-fb1833f50b87" +version = "0.2.1" + [[TableTraits]] deps = ["IteratorInterfaceExtensions"] git-tree-sha1 = "b1ad568ba658d8cbb3b892ed5380a6f3e781a81e" @@ -1048,10 +1179,10 @@ uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" version = "1.3.0" [[VectorizationBase]] -deps = ["CpuId", "LLVM", "Libdl", "LinearAlgebra"] -git-tree-sha1 = "bb72c58beab6c9e544851f5373fcd72f8f1f157a" +deps = ["CpuId", "Libdl", "LinearAlgebra"] +git-tree-sha1 = "1f2c3b362447e9dc031d09bf5a887b2ae01fcbd8" uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" -version = "0.12.21" +version = "0.12.23" [[VersionParsing]] git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f" @@ -1064,12 +1195,30 @@ git-tree-sha1 = "b9b450c99a3ca1cc1c6836f560d8d887bcbe356e" uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" version = "0.1.2" +[[Widgets]] +deps = ["Colors", "Dates", "Observables", "OrderedCollections"] +git-tree-sha1 = "fc0feda91b3fef7fe6948ee09bb628f882b49ca4" +uuid = "cc8bc4a8-27d6-5769-a93b-9d913e69aa62" +version = "0.6.2" + +[[WoodburyMatrices]] +deps = ["LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "28ffe06d28b1ba8fdb2f36ec7bb079fac81bac0d" +uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6" +version = "0.5.2" + [[Zlib_jll]] deps = ["Libdl", "Pkg"] git-tree-sha1 = "622d8b6dc0c7e8029f17127703de9819134d1b71" uuid = "83775a58-1f1d-513f-b197-d71354ab007a" version = "1.2.11+14" +[[Zygote]] +deps = ["DiffRules", "FFTW", "FillArrays", "ForwardDiff", "IRTools", "InteractiveUtils", "LinearAlgebra", "MacroTools", "NNlib", "NaNMath", "Random", "Requires", "SpecialFunctions", "Statistics", "ZygoteRules"] +git-tree-sha1 = "f8329b595c465caf3ca87c4f744e6041a4983e43" +uuid = "e88e6eb3-aa80-5325-afca-941959d7151f" +version = "0.4.8" + [[ZygoteRules]] deps = ["MacroTools"] git-tree-sha1 = "b3b4882cc9accf6731a08cc39543fbc6b669dca8" diff --git a/Project.toml b/Project.toml index 4d3a771d0..336bbfc15 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,9 @@ Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" @@ -17,5 +19,6 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/examples/GPEmulator/learn_noise.jl b/examples/GPEmulator/learn_noise.jl index 39d2b82c2..3796bb1c0 100644 --- a/examples/GPEmulator/learn_noise.jl +++ b/examples/GPEmulator/learn_noise.jl @@ -72,7 +72,6 @@ noise_samples = rand(MvNormal(μ, Σ), n)' # y = G(x) + η Y = gx .+ noise_samples -input_cov = cov(noise_samples, dims=1) # Fit 2D Gaussian Process regression model # (To be precise: We fit two models, one that predicts y1 from x1 and x2, and @@ -84,7 +83,7 @@ input_cov = cov(noise_samples, dims=1) # in the training phase via an optimization procedure. # Because of the SVD transformation applied to the output, we expect the # learned noise to be close to 1. -gpobj1 = GPObj(X, Y, input_cov, gppackage, GPkernel=nothing, +gpobj1 = GPObj(X, Y, gppackage, GPkernel=nothing, obs_noise_cov=Σ, normalized=true, noise_learn=true, prediction_type=pred_type) println("\n-----------") diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl index 0752f4e7a..7d2050f12 100644 --- a/examples/GPEmulator/plot_GP.jl +++ b/examples/GPEmulator/plot_GP.jl @@ -78,7 +78,6 @@ noise_samples = rand(MvNormal(μ, Σ), n)' # y = G(x) + η Y = gx .+ noise_samples -data_cov = cov(noise_samples, dims=1) # Fit 2D Gaussian Process regression model # (To be precise: We fit two models, one that predicts y1 from x1 and x2, @@ -87,7 +86,7 @@ data_cov = cov(noise_samples, dims=1) # exponential kernel. # Setting noise_learn=true leads to the addition of white noise to the # kernel -gpobj = GPObj(X, Y, data_cov, gppackage, GPkernel=nothing, +gpobj = GPObj(X, Y, gppackage, GPkernel=nothing, obs_noise_cov=Σ, normalized=true, noise_learn=true, prediction_type=pred_type) # Plot mean and variance of the predicted observables y1 and y2 diff --git a/src/EKI.jl b/src/EKI.jl index 455f10ceb..94b551935 100644 --- a/src/EKI.jl +++ b/src/EKI.jl @@ -44,7 +44,7 @@ end function EKIObj(parameters::Array{FT, 2}, parameter_names::Vector{String}, t_mean, - t_cov::Array{FT, 2}; + obs_noise_cov::Array{FT, 2}; Δt=FT(1)) where {FT<:AbstractFloat} # ensemble size @@ -60,7 +60,7 @@ function EKIObj(parameters::Array{FT, 2}, # timestep store Δt = Array([Δt]) - EKIObj{FT,IT}(u, parameter_names, t_mean, t_cov, N_ens, g, err, Δt) + EKIObj{FT,IT}(u, parameter_names, t_mean, obs_noise_cov, N_ens, g, err, Δt) end diff --git a/src/GPEmulator.jl b/src/GPEmulator.jl index 7ecf610aa..5fe10913c 100644 --- a/src/GPEmulator.jl +++ b/src/GPEmulator.jl @@ -17,6 +17,7 @@ export predict export GPJL, SKLJL export YType, FType +export svd_transform, svd_reverse_transform_mean_cov using Plots @@ -60,44 +61,46 @@ struct GPObj{FT<:AbstractFloat, GPM} data::Array{FT} "mean of input; 1 x input_dim" input_mean::Array{FT} - "square root of the inverse of the input covariance matrix; input_dim x input_dim" - sqrt_inv_input_cov::Union{Nothing, Array{FT, 2}} "the Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs" models::Vector + "the singular value decomposition of obs_noise_cov, such that obs_noise_cov = decomposition.U * Diagonal(decomposition.S) * decomposition.Vt." + decomposition::Union{SVD, Nothing} "whether to fit GP models on normalized inputs ((inputs - input_mean) * sqrt_inv_input_cov)" normalized::Bool - "the singular value decomposition matrix" - decomposition::SVD - "the normalization matrix in the svd transformed space" - sqrt_singular_values::Array{FT,2} + "square root of the inverse of the input covariance matrix; input_dim x input_dim" + sqrt_inv_input_cov::Union{Nothing, Array{FT, 2}} "prediction type (`y` to predict the data, `f` to predict the latent function)" prediction_type::Union{Nothing, PredictionType} end """ - GPObj(inputs, data, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true, svdflag=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} +GPObj(inputs::Array{FT, 2}, data::Array{FT, 2}, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, obs_noise_cov::Union{Array{FT, 2}, Nothing}=nothing, normalized::Bool=true, noise_learn::Bool=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} Inputs and data of size N_samples x input_dim (both arrays will be transposed in the construction of the GPObj) - - `GPkernel` - GaussianProcesses kernel object. If not supplied, a default - kernel is used. The default kernel is the sum of a Squared - Exponential kernel and white noise. + - `GPkernel` - GaussianProcesses kernel object. If not supplied, a default Squared Exponential kernel is used. """ function GPObj( inputs::Array{FT, 2}, data::Array{FT, 2}, - truth_cov::Array{FT, 2}, package::GPJL; GPkernel::Union{K, KPy, Nothing}=nothing, + obs_noise_cov::Union{Array{FT, 2}, Nothing}=nothing, normalized::Bool=true, noise_learn::Bool=true, prediction_type::PredictionType=YType()) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject} - # Consistency check - err = "Number of inputs is not equal to the number of data." - size(inputs, 1) == size(data, 1) || throw(ArgumentError(err)) + # Consistency checks + input_dim = size(inputs, 2) + output_dim = size(data, 2) + err1 = "Number of inputs is not equal to the number of data." + size(inputs, 1) == size(data, 1) || throw(ArgumentError(err1)) + if obs_noise_cov != nothing + err2 = "obs_noise_cov must be of size ($output_dim, $output_dim), got $(size(obs_noise_cov))" + size(obs_noise_cov) == (output_dim, output_dim) || throw(ArgumentError(err2)) + end # Initialize vector of GP models models = Any[] @@ -114,9 +117,6 @@ function GPObj( if normalized # Normalize and transpose (since the inputs have to be of size # input_dim x N_samples to pass to GPE()) - if ndims(inputs) == 1 - error("Cov not defined for 1d input; can't normalize") - end sqrt_inv_input_cov = convert(Array{FT}, sqrt(inv(cov(inputs, dims=1)))) GPinputs = convert(Array, ((inputs.-input_mean) * sqrt_inv_input_cov)') else @@ -124,16 +124,9 @@ function GPObj( GPinputs = convert(Array, inputs') end - # Transform the outputs (transformed_data will come out transposed from - # this transformation, i.e., with size output_dim x N_samples. This is also - # the size that will be required when passing the data to GPE(), so we leave - # transformed_data in that size. - # Create an SVD decomposition of the covariance: - # NB: svdfact() may be more efficient / provide positive definiteness information - # stored as: svd.U * svd.S *svd.Vt - decomposition = svd(truth_cov) - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - transformed_data = sqrt_singular_values_inv * decomposition.Vt * data' + # Transform data if obs_noise_cov available (if obs_noise_cov==nothing, + # transformed_data is equal to data) + transformed_data, decomposition = svd_transform(data, obs_noise_cov) # Use a default kernel unless a kernel was supplied to GPObj if GPkernel==nothing @@ -175,7 +168,7 @@ function GPObj( println("kernel in GPObj") GPkernel_i = deepcopy(kern) println(GPkernel_i) - GPdata_i = transformed_data[i, :] + GPdata_i = transformed_data[:, i] # GPE() arguments: # GPinputs: (input_dim x N_samples) # GPdata_i: (N_samples,) @@ -202,37 +195,41 @@ function GPObj( return GPObj{FT, typeof(package)}(inputs, data, input_mean, - sqrt_inv_input_cov, models, - normalized, decomposition, - inv(sqrt_singular_values_inv), + normalized, + sqrt_inv_input_cov, prediction_type) end """ - GPObj(inputs, data, package::SKLJL, GPkernel::Union{K, KPy, Nothing}=nothing, normalized::Bool=true) where {K<:Kernel, KPy<:PyObject} +GPObj(inputs::Array{FT, 2}, data::Array{FT, 2}, package::SKLJL; GPkernel::Union{K, KPy, Nothing}=nothing, obs_noise_cov::Union{Array{FT, 2}, Nothing}=nothing, normalized::Bool=true, noise_learn::Bool=true, prediction_type::PredictionType=YType()) where {K<:Kernel, KPy<:PyObject} -Inputs and data of size N_samples x input_dim (both arrays will be transposed in the construction of the GPObj) - - `GPkernel` - GaussianProcesses or ScikitLearn kernel object. If not supplied, - a default kernel is used. The default kernel is the sum of a - Squared Exponential kernel and white noise. +Inputs and data of size N_samples x input_dim + + - `GPkernel` - ScikitLearn kernel object. If not supplied, a default Squared Exponential kernel is used. """ function GPObj( inputs::Array{FT, 2}, data::Array{FT, 2}, - truth_cov::Array{FT, 2}, package::SKLJL; GPkernel::Union{K, KPy, Nothing}=nothing, + obs_noise_cov::Union{Array{FT, 2}, Nothing}=nothing, normalized::Bool=true, noise_learn::Bool=true, prediction_type::PredictionType=YType()) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject} - # Consistency check - err = "Number of inputs is not equal to the number of data." - size(inputs, 1) == size(data, 1) || throw(ArgumentError(err)) + # Consistency checks + input_dim = size(inputs, 2) + output_dim = size(data, 2) + err1 = "Number of inputs is not equal to the number of data." + size(inputs, 1) == size(data, 1) || throw(ArgumentError(err1)) + if obs_noise_cov != nothing + err2 = "obs_noise_cov must be of size ($output_dim, $output_dim), got $(size(obs_noise_cov))" + size(obs_noise_cov) == (output_dim, output_dim) || throw(ArgumentError(err2)) + end # Make sure inputs and data are arrays of type FT inputs = convert(Array{FT}, inputs) @@ -250,23 +247,15 @@ function GPObj( input_mean = reshape(mean(inputs, dims=1), 1, :) sqrt_inv_input_cov = nothing if normalized - if ndims(inputs) == 1 - error("Cov not defined for 1d input; can't normalize") - end sqrt_inv_input_cov = convert(Array{FT}, sqrt(inv(cov(inputs, dims=1)))) GPinputs = (inputs .- input_mean) * sqrt_inv_input_cov else GPinputs = inputs end - # Transform the outputs - # Create an SVD decomposition of the covariance: - # NB: svdfact() may be more efficient / provide positive definiteness information - # stored as: svd.U * svd.S *svd.Vt - decomposition = svd(truth_cov) - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - transformed_data = sqrt_singular_values_inv * decomposition.Vt * data' - + # Transform data if obs_noise_cov available (if obs_noise_cov==nothing, + # transformed_data is equal to data) + transformed_data, decomposition = svd_transform(data, obs_noise_cov) if GPkernel==nothing println("Using default squared exponential kernel, learning length scale and variance parameters") # Create default squared exponential kernel @@ -279,11 +268,11 @@ function GPObj( println("Using default squared exponential kernel:", kern) else kern = deepcopy(GPkernel) - println("Using user-defined kernel",kern) + println("Using user-defined kernel", kern) end if noise_learn - # Add white noise to kernel + # Add white noise to kernel white_noise_level = 1.0 white = WhiteKernel(noise_level=white_noise_level, noise_level_bounds=(1e-05, 10.0)) @@ -302,12 +291,11 @@ function GPObj( for i in 1:N_models GPkernel_i = deepcopy(kern) - GPdata_i = transformed_data[i, :] + GPdata_i = transformed_data[:, i] m = GaussianProcessRegressor(kernel=GPkernel_i, n_restarts_optimizer=10, alpha=regularization_noise) - println(m.kernel_) # ScikitLearn.fit! arguments: # GPinputs: (N_samples x input_dim) # GPdata_i: (N_samples,) @@ -316,41 +304,45 @@ function GPObj( println(m.kernel.hyperparameters) print("Completed training of: ") end - print(i,", ") + println(i, ", ") push!(models, m) - println(m.kernel_) + println(m.kernel) end return GPObj{FT, typeof(package)}(inputs, data, input_mean, - sqrt_inv_input_cov, models, - normalized, decomposition, - inv(sqrt_singular_values_inv), + normalized, + sqrt_inv_input_cov, YType()) end - """ - predict(gp::GPObj{FT,GPJL}, new_inputs::Array{FT, 2}) where {FT} = predict(gp, new_inputs, gp.prediction_type) + predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}, transform_to_real::Bool=false) where {FT} = predict(gp, new_inputs, gp.prediction_type) Evaluate the GP model(s) at new inputs. - `gp` - a GPObj - `new_inputs` - inputs for which GP model(s) is/are evaluated; N_new_inputs x input_dim +Returns the predicted mean(s) and covariance(s) at the input points. +Means: matrix of size N_new_inputs x output_dim +Covariances: vector of length N_new_inputs, each element is a matrix of size output_dim x output_dim. If the output is 1-dimensional, a N_new_inputs x 1 array of scalar variances is returned rather than a vector of 1x1 matrices. + Note: If gp.normalized == true, the new inputs are normalized prior to the prediction """ -predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}; transform_to_real=false) where {FT} = predict(gp, new_inputs, transform_to_real, gp.prediction_type) +predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}; transform_to_real::Bool=false) where {FT} = predict(gp, new_inputs, transform_to_real, gp.prediction_type) function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}, transform_to_real, ::FType) where {FT} + # Check if the size of new_inputs is consistent with the GP model's input # dimension. # new_inputs should be of size N_new_inputs x input_dim, so input_dim # has to equal the input dimension of the GP model(s). input_dim = size(gp.inputs, 2) output_dim = size(gp.data, 2) + N_new_inputs = size(new_inputs, 1) err = "GP object and input observations do not have consistent dimensions" size(new_inputs, 2) == input_dim || throw(ArgumentError(err)) @@ -359,43 +351,46 @@ function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}, transform_to_rea end M = length(gp.models) - μσ2 = [predict_f(gp.models[i], new_inputs) for i in 1:M] - # Return mean(s) and variance(s) - μ = reshape(vcat(first.(μσ2)...), size(new_inputs, 1), output_dim) - σ2 = reshape(vcat(last.(μσ2)...), size(new_inputs, 1), output_dim) - - if transform_to_real - # Revert back to the original (correlated) coordinate system - # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP - transformed_μT= gp.decomposition.V * gp.sqrt_singular_values * μ' - # Get the means back into size N_new_inputs x output_dim - transformed_μ = convert(Array, transformed_μT') - transformed_σ2 = zeros(size(σ2)) - # Back transformation of covariance - for j in 1:size(σ2, 1) - σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[j, :]) * gp.sqrt_singular_values * gp.decomposition.Vt - # Extract diagonal elements from σ2_j and add them to σ2 as column - # vectors - # TODO: Should we have an optional input argument `full_cov` - if - # set to true, the full covariance matrix is returned instead of - # only variances (default: false)? - transformed_σ2[j, :] = diag(σ2_j) - end - return transformed_μ, transformed_σ2 + # Predicts columns of inputs so must be transposed to size + # input_dim x N_new_inputs + new_inputs_transp = convert(Array{FT}, new_inputs') + μσ2 = [predict_f(gp.models[i], new_inputs_transp) for i in 1:M] + # Return mean(s) and variance(s) + μ = reshape(vcat(first.(μσ2)...), N_new_inputs, output_dim) + σ2 = reshape(vcat(last.(μσ2)...), N_new_inputs, output_dim) + + if transform_to_real && gp.decomposition != nothing + μ_pred, σ2_pred = svd_reverse_transform_mean_cov(μ, σ2, + gp.decomposition) + elseif transform_to_real && gp.decomposition == nothing + err = """Need SVD decomposition to transform back to original space, + but GPObj.decomposition == nothing. + Try setting transform_to_real=false""" + throw(ArgumentError(err)) else - # No back transformation needed - return μ, σ2 + μ_pred = μ + # Convert to vector of (diagonal) matrices to match the format of + # σ2_pred when transform_to_real=true + σ2_pred = vec([Diagonal(σ2[j, :]) for j in 1:N_new_inputs]) end + + if output_dim == 1 + σ2_pred = reshape([σ2_pred[i][1] for i in 1:N_new_inputs], N_new_inputs, 1) + end + + return μ_pred, σ2_pred end function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}, transform_to_real, ::YType) where {FT} + # Check if the size of new_inputs is consistent with the GP model's input # dimension. # new_inputs should be of size N_new_inputs x input_dim, so input_dim # has to equal the input dimension of the GP model(s). input_dim = size(gp.inputs, 2) output_dim = size(gp.data, 2) + N_new_inputs = size(new_inputs, 1) err = "GP object and input observations do not have consistent dimensions" size(new_inputs, 2) == input_dim || throw(ArgumentError(err)) @@ -407,39 +402,35 @@ function predict(gp::GPObj{FT, GPJL}, new_inputs::Array{FT, 2}, transform_to_rea # Predicts columns of inputs so must be transposed to size # input_dim x N_new_inputs new_inputs_transp = convert(Array{FT}, new_inputs') - μσ2 = [predict_y(gp.models[i], new_inputs_transp) for i in 1:M] + μσ2 = [predict_f(gp.models[i], new_inputs_transp) for i in 1:M] # Return mean(s) and variance(s) - μ = reshape(vcat(first.(μσ2)...), size(new_inputs, 1), output_dim) - σ2 = reshape(vcat(last.(μσ2)...), size(new_inputs, 1), output_dim) - - if transform_to_real - # Revert back to the original (correlated) coordinate system - # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP - transformed_μT= gp.decomposition.V * gp.sqrt_singular_values * μ' - # Get the means back into size N_new_inputs x output_dim - transformed_μ = convert(Array, transformed_μT') - transformed_σ2 = zeros(size(σ2)) - # Back transformation of covariance - for j in 1:size(σ2, 1) - σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[j, :]) * gp.sqrt_singular_values * gp.decomposition.Vt - # Extract diagonal elements from σ2_j and add them to σ2 as column - # vectors - # TODO: Should we have an optional input argument `full_cov` - if - # set to true, the full covariance matrix is returned instead of - # only variances (default: false)? - transformed_σ2[j, :] = diag(σ2_j) - end - return transformed_μ, transformed_σ2 - + μ = reshape(vcat(first.(μσ2)...), N_new_inputs, output_dim) + σ2 = reshape(vcat(last.(μσ2)...), N_new_inputs, output_dim) + + if transform_to_real && gp.decomposition != nothing + μ_pred, σ2_pred = svd_reverse_transform_mean_cov(μ, σ2, + gp.decomposition) + elseif transform_to_real && gp.decomposition == nothing + err = """Need SVD decomposition to transform back to original space, + but GPObj.decomposition == nothing. + Try setting transform_to_real=false""" + throw(ArgumentError(err)) else - # No back transformation needed - return μ, σ2 + μ_pred = μ + # Convert to vector of (diagonal) matrices to match the format of + # σ2_pred when transform_to_real=true + σ2_pred = vec([Diagonal(σ2[j, :]) for j in 1:N_new_inputs]) + end + + if output_dim == 1 + σ2_pred = reshape([σ2_pred[i][1] for i in 1:N_new_inputs], N_new_inputs, 1) end + return μ_pred, σ2_pred end -function predict(gp::GPObj{FT, SKLJL}, new_inputs::Array{FT, 2}, transform_to_real) where {FT} +function predict(gp::GPObj{FT, SKLJL}, new_inputs::Array{FT, 2}, transform_to_real::Bool=false) where {FT} # Check if the size of new_inputs is consistent with the GP model's input # dimension. @@ -447,6 +438,7 @@ function predict(gp::GPObj{FT, SKLJL}, new_inputs::Array{FT, 2}, transform_to_re # has to equal the input dimension of the GP model(s). input_dim = size(gp.inputs, 2) output_dim = size(gp.data, 2) + N_new_inputs = size(new_inputs, 1) err = "GP object and input observations do not have consistent dimensions" size(new_inputs, 2) == input_dim || throw(ArgumentError(err)) @@ -458,36 +450,107 @@ function predict(gp::GPObj{FT, SKLJL}, new_inputs::Array{FT, 2}, transform_to_re # Predicts rows of inputs; no need to transpose μσ = [gp.models[i].predict(new_inputs, return_std=true) for i in 1:M] # Return mean(s) and standard deviations(s) - μ = reshape(vcat(first.(μσ)...), size(new_inputs)[1], output_dim) - σ = reshape(vcat(last.(μσ)...), size(new_inputs)[1], output_dim) - σ2 = σ*σ - - if transform_to_real - # Revert back to the original (correlated) coordinate system - # We created meanvGP = Dinv*Vt*meanv so meanv = V*D*meanvGP - transformed_μT= gp.decomposition.V * gp.sqrt_singular_values * μ' - # Get the means back into size N_new_inputs x output_dim - transformed_μ = convert(Array, transformed_μT') - println("size transformed μ") - println(size(transformed_μ)) - transformed_σ2 = zeros(size(σ2)) - # Back transformation of covariance - for j in 1:size(σ2, 1) - σ2_j = gp.decomposition.V * gp.sqrt_singular_values * Diagonal(σ2[j, :]) * gp.sqrt_singular_values * gp.decomposition.Vt - # Extract diagonal elements from σ2_j and add them to σ2 as column - # vectors - # TODO: Should we have an optional input argument `full_cov` - if - # set to true, the full covariance matrix is returned instead of - # only variances (default: false)? - transformed_σ2[j, :] = diag(σ2_j) - end - return transformed_μ, transformed_σ2 + μ = reshape(vcat(first.(μσ)...), N_new_inputs, output_dim) + σ = reshape(vcat(last.(μσ)...), N_new_inputs, output_dim) + σ2 = σ .* σ + + if transform_to_real && gp.decomposition != nothing + μ_pred, σ2_pred = svd_reverse_transform_mean_cov(μ, σ2, + gp.decomposition) + elseif transform_to_real && gp.decomposition == nothing + err = """Need SVD decomposition to transform back to original space, + but GPObj.decomposition == nothing. + Try setting transform_to_real=false""" + throw(ArgumentError(err)) + else + μ_pred = μ + # Convert to vector of (diagonal) matrices to match the format of + # σ2_pred when transform_to_real=true + σ2_pred = vec([Diagonal(σ2[j, :]) for j in 1:N_new_inputs]) + end + + if output_dim == 1 + σ2_pred = reshape([σ2_pred[i][1] for i in 1:N_new_inputs], N_new_inputs, 1) + end + + return μ_pred, σ2_pred +end + +""" +svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, Nothing}) where {FT} + +Apply a singular value decomposition (SVD) to the data + - `data` - GP training data/targets; N_samples x output_dim + - `obs_noise_cov` - covariance of observational noise + +Returns the transformed data and the decomposition, which is a matrix +factorization of type LinearAlgebra.SVD. + +Note: If F::SVD is the factorization object, U, S, V and Vt can be obtained via +F.U, F.S, F.V and F.Vt, such that A = U * Diagonal(S) * Vt. The singular values +in S are sorted in descending order. +""" + +function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, Nothing}) where {FT} + if obs_noise_cov != nothing + decomposition = svd(obs_noise_cov) + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) + transformed_data_T = sqrt_singular_values_inv * decomposition.Vt * data' + transformed_data = convert(Array, transformed_data_T') else - # No back transformation needed - return μ, σ2 + decomposition = nothing + transformed_data = data end + return transformed_data, decomposition end +function svd_transform(data::Vector{FT}, obs_noise_cov::Union{Array{FT, 2}, Nothing}) where {FT} + if obs_noise_cov != nothing + decomposition = svd(obs_noise_cov) + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) + transformed_data = sqrt_singular_values_inv * decomposition.Vt * data + else + decomposition = nothing + transformed_data = data + end + + return transformed_data, decomposition end + + +""" +svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::{Array{FT, 2}, decomposition::SVD) where {FT} + +Transform the mean and covariance back to the original (correlated) coordinate system + - `μ` - predicted mean; N_predicted_points x output_dim + - `σ2` - predicted variance; N_predicted_points x output_dim + +Returns the transformed mean (N_predicted_points x output_dim) and variance. +Note that transforming the variance back to the original coordinate system +results in non-zero off-diagonal elements, so instead of just returning the +elements on the main diagonal (i.e., the variances), we return the full +covariance at each point, as a vector of length N_predicted_points, where +each element is a matrix of size output_dim x output_dim +""" + +function svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::Array{FT, 2}, decomposition::SVD) where {FT} + + N_predicted_points, output_dim = size(σ2) + # We created meanvGP = D_inv * Vt * mean_v so meanv = V * D * meanvGP + sqrt_singular_values= Diagonal(sqrt.(decomposition.S)) + transformed_μT = decomposition.V * sqrt_singular_values * μ' + # Get the means back into size N_prediction_points x output_dim + transformed_μ = convert(Array, transformed_μT') + transformed_σ2 = [zeros(output_dim, output_dim) for i in 1:N_predicted_points] + # Back transformation + for j in 1:N_predicted_points + σ2_j = decomposition.V * sqrt_singular_values * Diagonal(σ2[j, :]) * sqrt_singular_values * decomposition.Vt + transformed_σ2[j] = σ2_j + end + + return transformed_μ, transformed_σ2 +end + +end diff --git a/src/MCMC.jl b/src/MCMC.jl index 1d827f311..d9166daf3 100644 --- a/src/MCMC.jl +++ b/src/MCMC.jl @@ -32,9 +32,9 @@ struct MCMCObj{FT<:AbstractFloat, IT<:Int} "a single sample from the observations. Can e.g. be picked from an Obs struct using get_obs_sample" obs_sample::Vector{FT} "covariance of the observational noise" - obs_cov::Array{FT, 2} - "inverse of obs_cov" - obs_covinv::Array{FT, 2} + obs_noise_cov::Array{FT, 2} + "inverse of obs_noise_cov" + obs_noise_covinv::Array{FT, 2} "array of length N_parameters with the parameters' prior distributions" prior::Array{Prior, 1} "MCMC step size" @@ -57,7 +57,7 @@ end """ MCMCObj(obs_sample::Vector{FT}, - obs_cov::Array{FT, 2}, + obs_noise_cov::Array{FT, 2}, priors::Array{Prior, 1}, step::FT, param_init::Vector{FT}, @@ -70,7 +70,7 @@ where max_iter is the number of MCMC steps to perform (e.g., 100_000) """ function MCMCObj( obs_sample::Vector{FT}, - obs_cov::Array{FT, 2}, + obs_noise_cov::Array{FT, 2}, priors::Array{Prior, 1}, step::FT, param_init::Vector{FT}, @@ -85,9 +85,7 @@ function MCMCObj( # We need to transform obs_sample into the correct space if svdflag println("Applying SVD to decorrelating outputs, if not required set svdflag=false") - decomposition = svd(obs_cov) # svd.U * svd.S * svd.Vt (can also get V) - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) # diagonal matrix of 1/eigenvalues - obs_sample = sqrt_singular_values_inv * decomposition.Vt * obs_sample + obs_sample, unused = svd_transform(obs_sample, obs_noise_cov) else println("Assuming independent outputs.") end @@ -98,14 +96,14 @@ function MCMCObj( param = param_init_copy log_posterior = [nothing] iter = [1] - obs_covinv = inv(obs_cov) + obs_noise_covinv = inv(obs_noise_cov) accept = [0] if algtype != "rwm" error("only random walk metropolis 'rwm' is implemented so far") end MCMCObj{FT,IT}(obs_sample, - obs_cov, - obs_covinv, + obs_noise_cov, + obs_noise_covinv, priors, [step], burnin, @@ -169,7 +167,7 @@ function log_likelihood(mcmc::MCMCObj{FT}, log_rho = FT[0] if gvar == nothing diff = g - mcmc.obs_sample - log_rho[1] = -FT(0.5) * diff' * mcmc.obs_covinv * diff + log_rho[1] = -FT(0.5) * diff' * mcmc.obs_noise_covinv * diff else gcov_inv = inv(Diagonal(gvar)) log_gpfidelity = -FT(0.5) * log(det(Diagonal(gvar))) # = -0.5 * sum(log.(gvar)) @@ -226,10 +224,13 @@ function find_mcmc_step!(mcmc_test::MCMCObj{FT}, gpobj::GPObj{FT}) where {FT} local acc_ratio while mcmc_accept == false - param = convert(Array{FT, 2}, mcmc_test.param') + param = reshape(mcmc_test.param, 1, :) gp_pred, gp_predvar = predict(gpobj, param) - - mcmc_sample!(mcmc_test, vec(gp_pred), vec(gp_predvar)) + if ndims(gp_predvar[1]) != 0 + mcmc_sample!(mcmc_test, vec(gp_pred), diag(gp_predvar[1])) + else + mcmc_sample!(mcmc_test, vec(gp_pred), vec(gp_predvar)) + end it += 1 if it % 2000 == 0 countmcmc += 1 @@ -276,13 +277,15 @@ function sample_posterior!(mcmc::MCMCObj{FT,IT}, max_iter::IT) where {FT,IT<:Int} for mcmcit in 1:max_iter - param = convert(Array{FT, 2}, mcmc.param') - # test predictions (param' is 1 x N_parameters) + param = reshape(mcmc.param, 1, :) + # test predictions (param is 1 x N_parameters) gp_pred, gp_predvar = predict(gpobj, param) - gp_pred = cat(gp_pred..., dims=2) - gp_predvar = cat(gp_predvar..., dims=2) - mcmc_sample!(mcmc, vec(gp_pred), vec(gp_predvar)) + if ndims(gp_predvar[1]) != 0 + mcmc_sample!(mcmc, vec(gp_pred), diag(gp_predvar[1])) + else + mcmc_sample!(mcmc, vec(gp_pred), vec(gp_predvar)) + end end end diff --git a/src/Observations.jl b/src/Observations.jl index c1ba2b2ed..9d2dd2505 100644 --- a/src/Observations.jl +++ b/src/Observations.jl @@ -15,83 +15,140 @@ Structure that contains the observations $(DocStringExtensions.FIELDS) """ struct Obs{FT<:AbstractFloat} - "vector of observational samples, each of length N_data" + "vector of observational samples, each of length sample_dim" samples::Vector{Vector{FT}} - "covariance of the observational noise (assumed to be normally - distributed); N_data x N_data. If not supplied, cov is set to a - diagonal matrix whose non-zero elements are the sample variances" - cov::Union{Nothing, Array{FT, 2}} + "covariance of the observational noise (assumed to be normally + distributed); sample_dim x sample_dim (where sample_dim is the number of + elements in each sample), or a scalar if the sample dim is 1. If not + supplied, obs_noise_cov is set to a diagonal matrix whose non-zero elements + are the variances of the samples, or to a scalar variance in the case of + 1d samples. obs_noise_cov is set to nothing if only a single sample is + provided." + obs_noise_cov::Union{Array{FT, 2}, FT, Nothing} "sample mean" - mean::Vector{FT} + mean::Union{Vector{FT}, FT} "names of the data" - data_names::Vector{String} + data_names::Union{Vector{String}, String} end # Constructors function Obs(samples::Vector{Vector{FT}}, - data_names::Vector{String}) where {FT<:AbstractFloat} - N_samples = length(samples) + data_names::Union{Vector{String}, String}) where {FT<:AbstractFloat} - # convert to N_samples x N_data to determine sample covariance + N_samples = length(samples) + # convert to N_samples x sample_dim to determine sample covariance if N_samples == 1 # only one sample - this requires a bit more data massaging temp1 = convert(Array, reshape(hcat(samples...)', N_samples, :)) temp = dropdims(temp1, dims=1) samplemean = vec(mean(temp, dims=2)) - cov = nothing + obs_noise_cov = nothing else temp = convert(Array, reshape(hcat(samples...)', N_samples, :)) - samplemean = vec(mean(temp, dims=1)) - cov = Statistics.cov(temp .- samplemean) + if size(temp, 2) == 1 + # We have 1D samples, so the sample mean and covariance (which in + # this case is actually the covariance) are scalars + samplemean = mean(temp) + obs_noise_cov = var(temp) + else + samplemean = vec(mean(temp, dims=1)) + obs_noise_cov = cov(temp) + end end - Obs(samples, cov, samplemean, data_names) + Obs(samples, obs_noise_cov, samplemean, data_names) end -function Obs(samples::Array{FT}, - data_names::Vector{String}) where {FT<:AbstractFloat} - # samples is of size N_samples x N_data - if ndims(samples) == 1 +function Obs(samples::Array{FT, 2}, + data_names::Union{Vector{String}, String}) where {FT<:AbstractFloat} + + # samples is of size N_samples x sample_dim + N_samples, sample_dim = size(samples) + if N_samples == 1 # Only one sample, so there is no covariance to be computed and the # sample mean equals the sample itself - cov = nothing + obs_noise_cov = nothing samplemean = vec(samples) - samples = vec([vec(samples)]) + samples_vec = vec([vec(samples)]) else # convert matrix of samples to a vector of vectors - N_samples = size(samples, 1) - samplemean = vec(mean(samples, dims=1)) - cov = Statistics.cov(samples .- samplemean) - samples = [samples[i, :] for i in 1:N_samples] + samples_vec = [samples[i, :] for i in 1:N_samples] + if sample_dim == 1 + # We have 1D samples, so the sample mean and covariance (which in + # this case is actually the variance) are scalars + samplemean = mean(samples) + obs_noise_cov = var(samples) + else + samplemean = vec(mean(samples, dims=1)) + obs_noise_cov = cov(samples) + end end - Obs(samples, cov, samplemean, data_names) + Obs(samples_vec, obs_noise_cov, samplemean, data_names) end function Obs(samples::Vector{Vector{FT}}, - cov::Array{FT, 2}, - data_names::Vector{String}) where {FT<:AbstractFloat} + obs_noise_cov::Array{FT, 2}, + data_names::Union{Vector{String}, String}) where {FT<:AbstractFloat} + N_samples = length(samples) - # convert to N_samples x N_data to determine sample covariance + sample_dim = length(samples[1]) + # convert to N_samples x sample_dim to determine sample covariance if N_samples == 1 # only one sample - this requires a bit more data massaging temp1 = convert(Array, reshape(hcat(samples...)', N_samples, :)) temp = dropdims(temp1, dims=1) samplemean = mean(temp, dims=2) else - temp = convert(Array, reshape(hcat(samples...)', N_samples, :)) - samplemean = mean(temp, dims=1) + if sample_dim == 1 + # We have 1D samples, so the sample mean and covariance (which in + # this case is actually the covariance) are scalars + samplemean = mean(temp) + err = ("When sample_dim is 1, obs_cov_noise must be a scalar. + \tsample_dim: number of elements per observation sample") + @assert(ndims(obs_noise_cov) == 0, err) + else + temp = convert(Array, reshape(hcat(samples...)', N_samples, :)) + samplemean = vec(mean(temp, dims=1)) + err = ("obs_cov_noise must be of size sample_dim x sample_dim. + \tsample_dim: number of elements per observation sample") + @assert(size(obs_noise_cov) == (sample_dim, sample_dim), err) + end end - Obs(samples, cov, samplemean, data_names) + + + Obs(samples, obs_noise_cov, samplemean, data_names) end function Obs(samples::Array{FT, 2}, - cov::Array{FT, 2}, - data_names::Vector{String}) where {FT<:AbstractFloat} - # samples is of size N_samples x N_data - N_samples = size(samples, 1) - samplemean = vec(mean(samples, dims=1)) - # convert matrix of samples to a vector of vectors - samples = [samples[i, :] for i in 1:N_samples] - Obs(samples, cov, samplemean, data_names) + obs_noise_cov::Union{Array{FT, 2}, Nothing}, + data_names::Union{Vector{String}, String})where {FT<:AbstractFloat} + + # samples is of size N_samples x sample_dim + N_samples, sample_dim = size(samples) + if N_samples == 1 + # Only one sample, so there is no covariance to be computed and the + # sample mean equals the sample itself + obs_noise_cov = nothing + samplemean = vec(samples) + samples_vec = vec([vec(samples)]) + else + # convert matrix of samples to a vector of vectors + samples_vec = [samples[i, :] for i in 1:N_samples] + if sample_dim == 1 + # We have 1D samples, so the sample mean and covariance (which in + # this case is actually the variance) are scalars + samplemean = mean(samples) + err = ("When sample_dim is 1, obs_cov_noise must be a scalar. + \tsample_dim: number of elements per observation sample") + @assert(ndims(obs_noise_cov) == 0, err) + else + samplemean = vec(mean(samples, dims=1)) + err = ("obs_cov_noise must be of size sample_dim x sample_dim. + \tsample_dim: number of elements per observation sample") + @assert(size(obs_noise_cov) == (sample_dim, sample_dim), err) + end + end + + Obs(samples_vec, obs_noise_cov, samplemean, data_names) end end # module Observations diff --git a/test/Cloudy/runtests.jl b/test/Cloudy/runtests.jl index 3d218da2e..5b2eb27c0 100644 --- a/test/Cloudy/runtests.jl +++ b/test/Cloudy/runtests.jl @@ -10,7 +10,7 @@ using LinearAlgebra using StatsPlots using GaussianProcesses using Plots -using StatsPlots +using Random # Import Calibrate-Emulate-Sample modules using CalibrateEmulateSample.EKI @@ -21,6 +21,8 @@ using CalibrateEmulateSample.GModel using CalibrateEmulateSample.Utilities using CalibrateEmulateSample.Priors +rng_seed = 41 +Random.seed!(rng_seed) ### ### Define the (true) parameters and their priors @@ -117,7 +119,7 @@ N_ens = 50 # number of ensemble members N_iter = 5 # number of EKI iterations # initial parameters: N_ens x N_params initial_params = EKI.construct_initial_ensemble(N_ens, priors; rng_seed=6) -ekiobj = EKI.EKIObj(initial_params, param_names, truth.mean, truth.cov, Δt=1.0) +ekiobj = EKI.EKIObj(initial_params, param_names, truth.mean, truth.obs_noise_cov, Δt=1.0) # Initialize a ParticleDistribution with dummy parameters. The parameters # will then be set in run_G_ensemble @@ -167,9 +169,9 @@ GPkernel = kern1 + kern2 + white # Get training points u_tp, g_tp = Utilities.extract_GP_tp(ekiobj, N_iter) normalized = true -gpobj = GPEmulator.GPObj(u_tp, g_tp, Γy, gppackage; GPkernel=GPkernel, - normalized=normalized, noise_learn=false, - prediction_type=pred_type) +gpobj = GPEmulator.GPObj(u_tp, g_tp, gppackage; GPkernel=GPkernel, + obs_noise_cov=Γy, normalized=normalized, + noise_learn=false, prediction_type=pred_type) # Check how well the Gaussian Process regression predicts on the # true parameters diff --git a/test/EKI/runtests.jl b/test/EKI/runtests.jl index c9b92224f..552ed07b3 100644 --- a/test/EKI/runtests.jl +++ b/test/EKI/runtests.jl @@ -34,7 +34,7 @@ using CalibrateEmulateSample.Priors Γy = noise_level^2 * convert(Array, Diagonal(yt)) μ = zeros(length(yt)) for i in 1:n_samples - truth_samples[i, :] = yt + noise_level^2 * Distributions.rand(MvNormal(μ, Γy)) + truth_samples[i, :] = yt .+ Distributions.rand(MvNormal(μ, Γy)) end ### Calibrate: Ensemble Kalman Inversion diff --git a/test/GPEmulator/runtests.jl b/test/GPEmulator/runtests.jl index 143f76466..20873c043 100644 --- a/test/GPEmulator/runtests.jl +++ b/test/GPEmulator/runtests.jl @@ -1,7 +1,11 @@ # Import modules using Random -using GaussianProcesses using Test +using GaussianProcesses +using ScikitLearn +using Statistics +@sk_import gaussian_process : GaussianProcessRegressor +@sk_import gaussian_process.kernels : (RBF, WhiteKernel, ConstantKernel) using CalibrateEmulateSample.GPEmulator @@ -11,36 +15,126 @@ using CalibrateEmulateSample.GPEmulator rng_seed = 41 Random.seed!(rng_seed) + # ------------------------------------------------------------------------- + # Test case 1: 1D input, 1D output + # ------------------------------------------------------------------------- + # Training data - n = 10 # number of training points - x = 2.0 * π * rand(n) # predictors/features - y = sin.(x) + 0.05 * randn(n) # predictands/targets + n = 20 # number of training points + x = reshape(2.0 * π * rand(n), n, 1) # predictors/features: n x 1 + y = reshape(sin.(x) + 0.05 * randn(n), n, 1) # predictands/targets: n x 1 + + # Construct kernel: + # Squared exponential kernel (note that hyperparameters are on log scale) + # with observational noise + GPkernel = SE(log(1.0), log(1.0)) + + # These will be the test inputs at which predictions are made + new_inputs = reshape([0.0, π/2, π, 3*π/2, 2*π], 5, 1) + + # Fit Gaussian Process Regression models. + # GaussianProcesses.jl (GPJL) provides two predict functions, predict_y + # (which predicts the random variable y(θ)) and predict_y (which predicts + # the latent random variable f(θ)). + # ScikitLearn's Gaussian process regression (SKLJL) only offers one + # predict function, which predicts y. + # GPObj 1: GPJL, predict_y gppackage = GPJL() pred_type = YType() - # Construct kernel: - # Squared exponential kernel (note that hyperparameters are on log scale) - kern = SE(0.0, 0.0) - # log standard deviation of observational noise - logObsNoise = -1.0 - white = Noise(logObsNoise) - GPkernel = kern + white - - # Fit Gaussian Process Regression model - gpobj = GPObj(x, y, input_cov, gppackage; GPkernel=GPkernel, normalized=true, - noise_learn=true, prediction_type=pred_type) - new_inputs = [0.0, π/2, π, 3*π/2, 2*π] - μ, σ² = GPEmulator.predict(gpobj, new_inputs) - - @test_throws Exception GPObj(x, y, gppackage; GPkernel=GPkernel, - normalized=true, prediction_type=pred_type) - @test gpobj.inputs == x - @test gpobj.data == y - @test gpobj.input_mean[1] ≈ 3.524 atol=1e-2 - @test gpobj.sqrt_inv_input_cov == nothing - @test gpobj.prediction_type == pred_type - @test μ ≈ [0.0, 1.0, 0.0, -1.0, 0.0] atol=1.0 - @test σ² ≈ [0.859, 0.591, 0.686, 0.645, 0.647] atol=1e-2 + gpobj1 = GPObj(x, y, gppackage; GPkernel=GPkernel, obs_noise_cov=nothing, + normalized=false, noise_learn=true, + prediction_type=pred_type) + μ1, σ1² = GPEmulator.predict(gpobj1, new_inputs) + + @test gpobj1.inputs == x + @test gpobj1.data == y + @test gpobj1.input_mean[1] ≈ 3.0 atol=1e-2 + @test gpobj1.sqrt_inv_input_cov == nothing + @test gpobj1.prediction_type == pred_type + @test vec(μ1) ≈ [0.0, 1.0, 0.0, -1.0, 0.0] atol=0.3 + @assert(size(μ1)==(5,1)) + @test vec(σ1²) ≈ [0.017, 0.003, 0.004, 0.004, 0.009] atol=1e-2 + + # Check if normalization works + gpobj1_norm = GPObj(x, y, gppackage; GPkernel=GPkernel, + obs_noise_cov=nothing, normalized=true, + noise_learn=true, prediction_type=pred_type) + @test gpobj1_norm.sqrt_inv_input_cov ≈ [sqrt(1.0 / Statistics.var(x))] atol=1e-4 + + + # GPObj 2: GPJL, predict_f + pred_type = FType() + + gpobj2 = GPObj(x, y, gppackage; GPkernel=GPkernel, obs_noise_cov=nothing, + normalized=false, noise_learn=true, + prediction_type=pred_type) + μ2, σ2² = GPEmulator.predict(gpobj2, new_inputs) + # predict_y and predict_f should give the same mean + @test μ2 ≈ μ1 atol=1e-6 + + + # GPObj 3: SKLJL + gppackage = SKLJL() + pred_type = YType() + var = ConstantKernel(constant_value=1.0) + se = RBF(1.0) + GPkernel = var * se + + gpobj3 = GPObj(x, y, gppackage; GPkernel=GPkernel, obs_noise_cov=nothing, + normalized=false, noise_learn=true, + prediction_type=pred_type) + μ3, σ3² = GPEmulator.predict(gpobj3, new_inputs) + @test vec(μ3) ≈ [0.0, 1.0, 0.0, -1.0, 0.0] atol=0.3 + @test vec(σ3²) ≈ [0.016, 0.002, 0.003, 0.004, 0.003] atol=1e-2 + + + # ------------------------------------------------------------------------- + # Test case 2: 2D input, 2D output + # ------------------------------------------------------------------------- + + gppackage = GPJL() + pred_type = YType() + + # Generate training data + p = 2 # input dim + d = 2 # output dim + X = 2.0 * π * rand(n, p) + + # G(x1, x2) + g1x = sin.(X[:, 1]) .+ cos.(X[:, 2]) + g2x = sin.(X[:, 1]) .- cos.(X[:, 2]) + gx = [g1x g2x] + + # Add noise η + μ = zeros(d) + Σ = 0.1 * [[0.8, 0.2] [0.2, 0.5]] # d x d + noise_samples = rand(MvNormal(μ, Σ), n)' + + # y = G(x) + η + Y = gx .+ noise_samples + + transformed_Y, decomposition = svd_transform(Y, Σ) + @test size(transformed_Y) == size(Y) + transformed_Y, decomposition = svd_transform(Y[1, :], Σ) + @test size(transformed_Y) == size(Y[1, :]) + + gpobj4 = GPObj(X, Y, gppackage, GPkernel=nothing, obs_noise_cov=Σ, + normalized=true, noise_learn=true, + prediction_type=pred_type) + + new_inputs = zeros(4, 2) + new_inputs[2, :] = [π/2, π] + new_inputs[3, :] = [π, π/2] + new_inputs[4, :] = [3*π/2, 2*π] + μ4, σ4² = GPEmulator.predict(gpobj4, new_inputs, transform_to_real=true) + @test μ4[1, :] ≈ [1.0, -1.0] atol=0.5 + @test μ4[2, :] ≈ [0.0, 2.0] atol=0.5 + @test μ4[3, :] ≈ [0.0, 0.0] atol=0.5 + @test μ4[4, :] ≈ [0.0, -2.0] atol=0.5 + @test length(σ4²) == size(new_inputs, 1) + @test size(σ4²[1]) == (d, d) + end diff --git a/test/MCMC/runtests.jl b/test/MCMC/runtests.jl index de408065c..f192615ce 100644 --- a/test/MCMC/runtests.jl +++ b/test/MCMC/runtests.jl @@ -16,32 +16,28 @@ using CalibrateEmulateSample.GPEmulator # We need a GPObj to run MCMC, so let's reconstruct the one that's tested # in test/GPEmulator/runtests.jl - # Training data - n = 10 # number of training points - x = 2.0 * π * rand(n) # predictors/features - y = sin.(x) + 0.05 * randn(n) # predictands/targets + n = 20 # number of training points + x = reshape(2.0 * π * rand(n), n, 1) # predictors/features: n x 1 + σ2_y = reshape([0.05], 1, 1) + y = reshape(sin.(x) + rand(Normal(0, σ2_y[1]), n), n, 1) # predictands/targets: n x 1 gppackage = GPJL() pred_type = YType() # Construct kernel: # Squared exponential kernel (note that hyperparameters are on log scale) - kern = SE(0.0, 0.0) - # log standard deviation of observational noise - logObsNoise = -1.0 - white = Noise(logObsNoise) - GPkernel = kern + white + # with observational noise + GPkernel = SE(log(1.0), log(1.0)) - # Fit Gaussian Process Regression model - gpobj = GPObj(x, y, gppackage; GPkernel=GPkernel, normalized=false, - prediction_type=pred_type) + gpobj = GPObj(x, y, gppackage; GPkernel=GPkernel, obs_noise_cov=σ2_y, + normalized=false, noise_learn=true, + prediction_type=pred_type) ### Define prior umin = -1.0 umax = 6.0 prior = [Priors.Prior(Uniform(umin, umax), "u")] # prior on u obs_sample = [1.0] - cov = convert(Array, Diagonal([1.0])) # MCMC parameters mcmc_alg = "rwm" # random walk Metropolis @@ -51,7 +47,7 @@ using CalibrateEmulateSample.GPEmulator burnin = 0 step = 0.5 # first guess max_iter = 5000 - mcmc_test = MCMCObj(obs_sample, cov, prior, step, param_init, max_iter, + mcmc_test = MCMCObj(obs_sample, σ2_y, prior, step, param_init, max_iter, mcmc_alg, burnin) new_step = find_mcmc_step!(mcmc_test, gpobj) @@ -60,20 +56,22 @@ using CalibrateEmulateSample.GPEmulator max_iter = 100_000 # Now begin the actual MCMC - mcmc = MCMCObj(obs_sample, cov, prior, step, param_init, max_iter, - mcmc_alg, burnin) + mcmc = MCMCObj(obs_sample, σ2_y, prior, step, param_init, max_iter, + mcmc_alg, burnin, svdflag=true) sample_posterior!(mcmc, gpobj, max_iter) posterior = get_posterior(mcmc) post_mean = mean(posterior, dims=1)[1] - - @test mcmc.obs_sample == obs_sample - @test mcmc.obs_cov == cov - @test mcmc.obs_covinv == inv(cov) + + # We had svdflag=true, so the MCMCObj stores a transformed sample, + # 1.0/sqrt(0.05) * obs_sample ≈ 4.472 + @test mcmc.obs_sample ≈ [4.472] atol=1e-2 + @test mcmc.obs_noise_cov == σ2_y + @test mcmc.obs_noise_covinv == inv(σ2_y) @test mcmc.burnin == burnin @test mcmc.algtype == mcmc_alg @test mcmc.iter[1] == max_iter + 1 @test size(posterior) == (max_iter - burnin + 1, length(param_init)) - @test_throws Exception MCMCObj(obs_sample, cov, prior, step, param_init, + @test_throws Exception MCMCObj(obs_sample, σ2_y, prior, step, param_init, max_iter, "gibbs", burnin) @test post_mean ≈ π/2 atol=3e-1 diff --git a/test/Observations/runtests.jl b/test/Observations/runtests.jl index 02a6daa70..c2da502ec 100644 --- a/test/Observations/runtests.jl +++ b/test/Observations/runtests.jl @@ -7,27 +7,47 @@ using CalibrateEmulateSample.Observations @testset "Observations" begin # Generate samples as vector of vectors - n_data = 3 - n_samples = 5 - samples = [vec(i*ones(n_data)) for i in 1:n_samples] - data_names = ["d$(string(i))" for i in 1:n_data] + sample_dim = 3 # sample dimension, i.e., number of elements per sample + n_samples = 5 # number of samples + samples = [vec(i*ones(sample_dim)) for i in 1:n_samples] + data_names = ["d$(string(i))" for i in 1:sample_dim] obs = Obs(samples, data_names) @test obs.data_names == data_names @test obs.mean == [3.0, 3.0, 3.0] - @test obs.cov == 2.5 * ones(3, 3) + @test obs.obs_noise_cov == 2.5 * ones(3, 3) - # Generate a single sample - sample = [1.0, 2.0, 3.0] - obs = Obs(sample, data_names) - @test obs.mean == sample - @test obs.cov == nothing + # Generate samples as vector of vectors, pass a covariance to Obs + covar = ones(sample_dim, sample_dim) + obs = Obs(samples, covar, data_names) + @test obs.obs_noise_cov == covar + covar_wrong_dims = ones(sample_dim, sample_dim-1) + @test_throws AssertionError Obs(samples, covar_wrong_dims, data_names) # Generate samples as a 2d-array (each row corresponds to 1 sample) - samples = vcat([i*ones(n_data)' for i in 1:n_samples]...) + samples = vcat([i*ones(sample_dim)' for i in 1:n_samples]...) obs = Obs(samples, data_names) @test obs.mean == [3.0, 3.0, 3.0] - @test obs.cov == 2.5 * ones(3, 3) + @test obs.obs_noise_cov == 2.5 * ones(3, 3) + + # Generate samples as a 2d-array (each row corresponds to 1 sample), + # pass a covariance to Obs + obs = Obs(samples, covar, data_names) + @test obs.obs_noise_cov == covar + @test_throws AssertionError Obs(samples, covar_wrong_dims, data_names) + # Generate a single sample (a row vector) + sample = reshape([1.0, 2.0, 3.0], 1, 3) + obs = Obs(sample, data_names) + @test obs.mean == vec(sample) + @test obs.obs_noise_cov == nothing + + # Generate 1D-samples (a column vector) -- this should result in scalar + # values for the mean and obs_nosie_cov + sample = reshape([1.0, 2.0, 3.0], 3, 1) + data_name = "d1" + obs = Obs(sample, data_name) + @test obs.mean == 2.0 + @test obs.obs_noise_cov == 1.0 end From ca449560cb9c3442f076b957522d2159ff8fe809 Mon Sep 17 00:00:00 2001 From: Melanie Bieli Date: Thu, 16 Jul 2020 10:26:47 -0400 Subject: [PATCH 15/18] Update examples to work with updated code. --- examples/GPEmulator/learn_noise.jl | 8 ++++---- examples/GPEmulator/plot_GP.jl | 14 ++++++++++++-- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/examples/GPEmulator/learn_noise.jl b/examples/GPEmulator/learn_noise.jl index 3796bb1c0..64e78a199 100644 --- a/examples/GPEmulator/learn_noise.jl +++ b/examples/GPEmulator/learn_noise.jl @@ -54,7 +54,7 @@ pred_type = YType() # The observables y are related to the parameters x by: # y = G(x1, x2) + η, # where G(x1, x2) := [sin(x1) + cos(x2), sin(x1) - cos(x2)], and η ~ N(0, Σ) -n = 50 # number of training points +n = 200 # number of training points p = 2 # input dim d = 2 # output dim @@ -96,7 +96,7 @@ println("Learned noise parameter, σ_1:") learned_σ_1 = sqrt(gpobj1.models[1].kernel.kright.σ2) println("σ_1 = $learned_σ_1") # Check if the learned noise is approximately 1 -@assert(isapprox(learned_σ_1, 1.0; atol=0.07)) +@assert(isapprox(learned_σ_1, 1.0; atol=0.1)) println("\nKernel of the GP trained to predict y2 from x=(x1, x2):") println(gpobj1.models[2].kernel) @@ -104,13 +104,13 @@ println("Learned noise parameter, σ_2:") learned_σ_2 = sqrt(gpobj1.models[2].kernel.kright.σ2) println("σ_2 = $learned_σ_2") # Check if the learned noise is approximately 1 -@assert(isapprox(learned_σ_2, 1.0; atol=0.07)) +@assert(isapprox(learned_σ_2, 1.0; atol=0.1)) println("------------------------------------------------------------------\n") # For comparison: When noise_learn is set to false, the observational noise # is set to 1.0 and is not learned/optimized during the training. But thanks # to the SVD, 1.0 is the correct value to use. -gpobj2 = GPObj(X, Y, input_cov, gppackage, GPkernel=nothing, +gpobj2 = GPObj(X, Y, gppackage, GPkernel=nothing, obs_noise_cov=Σ, normalized=true, noise_learn=false, prediction_type=pred_type) println("\n-----------") diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl index 7d2050f12..4d17ce62e 100644 --- a/examples/GPEmulator/plot_GP.jl +++ b/examples/GPEmulator/plot_GP.jl @@ -3,7 +3,7 @@ using Random using GaussianProcesses using Distributions using Statistics -using Plots; pyplot(size = (1500, 700)) +using Plots; pyplot(size=(1500, 700)) Plots.scalefontsizes(1.3) using LinearAlgebra using CalibrateEmulateSample.GPEmulator @@ -98,10 +98,20 @@ X1, X2 = meshgrid(x1, x2) # Input for predict has to be of size N_samples x input_dim inputs = hcat(X1[:], X2[:]) +font = Plots.font("Helvetica", 18) +fontdict = Dict(:guidefont=>font, :xtickfont=>font, :ytickfont=>font, + :legendfont=>font) for y_i in 1:d - gp_mean, gp_var = GPEmulator.predict(gpobj, + # Predict on the grid points (note that `predict` returns the full + # covariance matrices, not just the variance -- gp_cov is a vector + # of covariance matrices) + gp_mean, gp_cov = GPEmulator.predict(gpobj, inputs, transform_to_real=true) + # Reshape gp_cov to size N_samples x output_dim + gp_var_temp = [diag(gp_cov[j]) for j in 1:length(gp_cov)] # (40000,) + gp_var = vcat([x' for x in gp_var_temp]...) # 40000 x 2 + mean_grid = reshape(gp_mean[:, y_i], n_pts, n_pts) p1 = plot(x1, x2, mean_grid, st=:surface, camera=(-30, 30), c=:cividis, xlabel="x1", ylabel="x2", zlabel="mean of y"*string(y_i), From f337bd93e59d595142e75722b0c8d07cbb05a9ad Mon Sep 17 00:00:00 2001 From: Melanie Bieli Date: Fri, 17 Jul 2020 19:22:35 -0400 Subject: [PATCH 16/18] Add plot of the difference between the prediction and the truth to `plot_GP.jl` Also: removed the import of GaussianProcesses.jl (which is not used) and changed the covariance used in that example to a diagonal matrix, as it doesn't really make sense to plot the variance otherwise (i.e., for full covariance matrices with off-diagonal elements). --- examples/GPEmulator/plot_GP.jl | 45 +++++++++++++++++++++++++++------- 1 file changed, 36 insertions(+), 9 deletions(-) diff --git a/examples/GPEmulator/plot_GP.jl b/examples/GPEmulator/plot_GP.jl index 4d17ce62e..e4be4ed98 100644 --- a/examples/GPEmulator/plot_GP.jl +++ b/examples/GPEmulator/plot_GP.jl @@ -1,6 +1,5 @@ # Import modules using Random -using GaussianProcesses using Distributions using Statistics using Plots; pyplot(size=(1500, 700)) @@ -22,7 +21,7 @@ using CalibrateEmulateSample.GPEmulator # and one that predicts y_i[2] from x_i. Fitting separate models for # # y_i[1] and y_i[2] requires the outputs to be independent - since this # # cannot be assumed to be the case a priori, the training data are # -# decorrelated by perfoming an Singulat Value Decomposition (SVD) on the # +# decorrelated by perfoming a Singular Value Decomposition (SVD) on the # # noise covariance Σ, and each model is then trained in the decorrelated # # space. # # # @@ -73,7 +72,7 @@ gx = [g1x g2x] # Add noise η μ = zeros(d) -Σ = 0.1 * [[0.8, 0.2] [0.2, 0.5]] # d x d +Σ = 0.1 * [[0.8, 0.0] [0.0, 0.5]] # d x d noise_samples = rand(MvNormal(μ, Σ), n)' # y = G(x) + η @@ -90,7 +89,7 @@ gpobj = GPObj(X, Y, gppackage, GPkernel=nothing, obs_noise_cov=Σ, normalized=true, noise_learn=true, prediction_type=pred_type) # Plot mean and variance of the predicted observables y1 and y2 -# For this, we generate test points on a x1-x2 grid +# For this, we generate test points on a x1-x2 grid. n_pts = 200 x1 = range(0.0, stop=2*π, length=n_pts) x2 = range(0.0, stop=2*π, length=n_pts) @@ -126,18 +125,46 @@ for y_i in 1:d savefig("GP_test_y"*string(y_i)*"_predictions.png") end -# Make plots of the true components of G(x1, x2) +# Plot the true components of G(x1, x2) g1_true = sin.(inputs[:, 1]) .+ cos.(inputs[:, 2]) -g1_grid = reshape(g1_true, n_pts, n_pts) -p3 = plot(x1, x2, g1_grid, st=:surface, camera=(-30, 30), c=:cividis, +g1_true_grid = reshape(g1_true, n_pts, n_pts) +p3 = plot(x1, x2, g1_true_grid, st=:surface, camera=(-30, 30), c=:cividis, xlabel="x1", ylabel="x2", zlabel="sin(x1) + cos(x2)", zguidefontrotation=90) savefig("GP_test_true_g1.png") g2_true = sin.(inputs[:, 1]) .- cos.(inputs[:, 2]) -g2_grid = reshape(g2_true, n_pts, n_pts) -p4 = plot(x1, x2, g2_grid, st=:surface, camera=(-30, 30), c=:cividis, +g2_true_grid = reshape(g2_true, n_pts, n_pts) +p4 = plot(x1, x2, g2_true_grid, st=:surface, camera=(-30, 30), c=:cividis, xlabel="x1", ylabel="x2", zlabel="sin(x1) - cos(x2)", zguidefontrotation=90) +g_true_grids = [g1_true_grid, g2_true_grid] savefig("GP_test_true_g2.png") + + +# Plot the difference between the truth and the mean of the predictions +for y_i in 1:d + # Predict on the grid points (note that `predict` returns the full + # covariance matrices, not just the variance -- gp_cov is a vector + # of covariance matrices) + gp_mean, gp_cov = GPEmulator.predict(gpobj, + inputs, + transform_to_real=true) + # Reshape gp_cov to size N_samples x output_dim + gp_var_temp = [diag(gp_cov[j]) for j in 1:length(gp_cov)] # (40000,) + gp_var = vcat([x' for x in gp_var_temp]...) # 40000 x 2 + + mean_grid = reshape(gp_mean[:, y_i], n_pts, n_pts) + var_grid = reshape(gp_var[:, y_i], n_pts, n_pts) + # Compute and plot 1/variance * (truth - prediction)^2 + zlabel = "1/var * (true_y"*string(y_i)*" - predicted_y"*string(y_i)*")^2" + p5 = plot(x1, x2, 1.0 ./ var_grid .* (g_true_grids[y_i] .- mean_grid).^2, + st=:surface, camera=(-30, 30), c=:magma, zlabel=zlabel, + xlabel="x1", ylabel="x2", + zguidefontrotation=90) + + savefig("GP_test_y"*string(y_i)*"_difference_truth_prediction.png") +end + + Plots.scalefontsizes(1/1.3) From 24ed6def6de0f6ca7f2f46b9adbe7eb84f464f8c Mon Sep 17 00:00:00 2001 From: Melanie Bieli Date: Sun, 19 Jul 2020 21:01:14 -0400 Subject: [PATCH 17/18] Remove GaussianProcesses dependency from `learn_noise.jl` Removed GaussianProcesses from `plot_GP.jl` in the previous commit, but forgot to also do it for `learn_noise.jl`. --- examples/GPEmulator/learn_noise.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/GPEmulator/learn_noise.jl b/examples/GPEmulator/learn_noise.jl index 64e78a199..496754343 100644 --- a/examples/GPEmulator/learn_noise.jl +++ b/examples/GPEmulator/learn_noise.jl @@ -1,6 +1,5 @@ # Import modules using Random -using GaussianProcesses using Distributions using Statistics using LinearAlgebra From 442269ce1794cfda96e4ee8cb0f056a569680448 Mon Sep 17 00:00:00 2001 From: odunbar Date: Tue, 21 Jul 2020 14:53:27 +0100 Subject: [PATCH 18/18] resolve CI fail --- test/Observations/runtests.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/Observations/runtests.jl b/test/Observations/runtests.jl index c14e811f8..6817215cc 100644 --- a/test/Observations/runtests.jl +++ b/test/Observations/runtests.jl @@ -52,5 +52,3 @@ using CalibrateEmulateSample.Observations @test obs.obs_noise_cov == 1.0 end - -end