From 5d9def724ba9795e8a6d2cff9ccda6057943dfe7 Mon Sep 17 00:00:00 2001 From: Thomas Jackson Date: Wed, 26 Jan 2022 00:04:01 -0500 Subject: [PATCH 1/4] Substitute types with abstractions Implements abstract typing done in EnsembleKalmanProcesses.jl#100. Adds support & tests for covariance as a dense matrix, Diagonal or UniformScaling type. Use abstract types in Emulator Combine logic for GaussianProcess predict() Use abstract types in GaussianProcess Use abstract types in MarkovChainMonteCarlo Add UniformScaling as a covariance type Correct SVD transformation Qualify use of GaussianProcesses.GPE, for clarity --- docs/src/API/Emulators.md | 1 - examples/Lorenz/Lorenz_example.jl | 6 +- examples/deprecated/Cloudy/Cloudy_example.jl | 5 +- src/Emulator.jl | 252 ++++++++----------- src/GaussianProcess.jl | 99 +++----- src/MarkovChainMonteCarlo.jl | 147 ++++++----- test/Emulator/runtests.jl | 14 +- 7 files changed, 238 insertions(+), 286 deletions(-) diff --git a/docs/src/API/Emulators.md b/docs/src/API/Emulators.md index f4c97aa11..12fcf7b25 100644 --- a/docs/src/API/Emulators.md +++ b/docs/src/API/Emulators.md @@ -5,7 +5,6 @@ CurrentModule = CalibrateEmulateSample.Emulators ``` ```@docs -Decomposition Emulator optimize_hyperparameters!(::Emulator) predict diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl index d1fe4955f..f31f7699f 100644 --- a/examples/Lorenz/Lorenz_example.jl +++ b/examples/Lorenz/Lorenz_example.jl @@ -317,11 +317,9 @@ optimize_hyperparameters!(emulator) # true parameters #if truncate_svd==1.0 if log_normal==false - y_mean, y_var = Emulators.predict(emulator, reshape(params_true, :, 1), - transform_to_real=true) + y_mean, y_var = Emulators.predict(emulator, params_true, transform_to_real=true) else - y_mean, y_var = Emulators.predict(emulator, reshape(log.(params_true), :, 1), - transform_to_real=true) + y_mean, y_var = Emulators.predict(emulator, log.(params_true), transform_to_real=true) end println("GP prediction on true parameters: ") diff --git a/examples/deprecated/Cloudy/Cloudy_example.jl b/examples/deprecated/Cloudy/Cloudy_example.jl index 3febf7877..f4223540c 100644 --- a/examples/deprecated/Cloudy/Cloudy_example.jl +++ b/examples/deprecated/Cloudy/Cloudy_example.jl @@ -225,10 +225,7 @@ optimize_hyperparameters!(emulator) # Check how well the Gaussian Process regression predicts on the # true parameters -y_mean, y_var = Emulators.predict( - emulator, reshape(transformed_params_true, :, 1); - transform_to_real=true -) +y_mean, y_var = Emulators.predict(emulator, transformed_params_true; transform_to_real=true) println("GP prediction on true parameters: ") println(vec(y_mean)) println("true data: ") diff --git a/src/Emulator.jl b/src/Emulator.jl index 5c7ca9c03..d3fca75ea 100644 --- a/src/Emulator.jl +++ b/src/Emulator.jl @@ -7,30 +7,10 @@ using LinearAlgebra using DocStringExtensions export Emulator -export Decomposition export optimize_hyperparameters! export predict -# SVD decomposition structure -""" - Decomposition{<:AbstractFloat, <:Int} - -Struct of SVD decomposition, containing (V,S,Vt), and the size of S, N. -""" -struct Decomposition{FT<:AbstractFloat, IT<:Int} - V::Array{FT,2} - Vt::Array{FT,2} - S::Array{FT} - N::IT -end - -# SVD decomposition constructor -function Decomposition(svd::SVD) - # svd.V is of type adjoint, transformed to Array with [:,:] - return Decomposition(svd.V[:,:], svd.Vt, svd.S, size(svd.S)[1]) -end - """ MachineLearningTool @@ -61,35 +41,35 @@ include("GaussianProcess.jl") #for GaussianProcess Structure used to represent a general emulator: """ -struct Emulator{FT} +struct Emulator{FT<:AbstractFloat} "Machine learning tool, defined as a struct of type MachineLearningTool" machine_learning_tool::MachineLearningTool - "normalized, standardized, transformed pairs given the Boolean's normalize_inputs, standardize_outputs, truncate_svd " + "normalized, standardized, transformed pairs given the Boolean's normalize_inputs, standardize_outputs, truncate_svd " training_pairs::PairedDataContainer{FT} - "mean of input; 1 × input_dim" - input_mean::Array{FT} + "mean of input; length input_dim" + input_mean::AbstractVector{FT} "square root of the inverse of the input covariance matrix; input_dim × input_dim" normalize_inputs::Bool "whether to fit models on normalized outputs outputs / standardize_outputs_factor" - sqrt_inv_input_cov::Union{Nothing, Array{FT, 2}} + sqrt_inv_input_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}, Nothing} " if normalizing: whether to fit models on normalized inputs ((inputs - input_mean) * sqrt_inv_input_cov)" standardize_outputs::Bool "if standardizing: Standardization factors (characteristic values of the problem)" - standardize_outputs_factors + standardize_outputs_factors::Union{AbstractVector{FT}, Nothing} "the singular value decomposition of obs_noise_cov, such that obs_noise_cov = decomposition.U * Diagonal(decomposition.S) * decomposition.Vt. NB: the svd may be reduced in dimensions" - decomposition::Union{Decomposition, Nothing} + decomposition::Union{SVD, Nothing} end # Constructor for the Emulator Object function Emulator( machine_learning_tool::MachineLearningTool, input_output_pairs::PairedDataContainer{FT}; - obs_noise_cov=nothing, + obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}, Nothing} = nothing, normalize_inputs::Bool = true, standardize_outputs::Bool = false, - standardize_outputs_factors::Union{Array{FT,1}, Nothing}=nothing, + standardize_outputs_factors::Union{AbstractVector{FT}, Nothing} = nothing, truncate_svd::FT=1.0 - ) where {FT<:AbstractFloat} +) where {FT<:AbstractFloat} # For Consistency checks input_dim, output_dim = size(input_output_pairs, 1) @@ -100,10 +80,10 @@ function Emulator( # [1.] Normalize the inputs? - input_mean = reshape(mean(get_inputs(input_output_pairs), dims=2), :, 1) #column vector + input_mean = vec(mean(get_inputs(input_output_pairs), dims=2)) #column vector sqrt_inv_input_cov = nothing if normalize_inputs - # Normalize (NB the inputs have to be of) size [input_dim × N_samples] to pass to GPE()) + # Normalize (NB the inputs have to be of) size [input_dim × N_samples] to pass to ML tool sqrt_inv_input_cov = sqrt(inv(Symmetric(cov(get_inputs(input_output_pairs), dims=2)))) training_inputs = normalize( get_inputs(input_output_pairs), @@ -143,14 +123,16 @@ function Emulator( # [4.] build an emulator build_models!(machine_learning_tool,training_pairs) - return Emulator{FT}(machine_learning_tool, - training_pairs, - input_mean, - normalize_inputs, - sqrt_inv_input_cov, - standardize_outputs, - standardize_outputs_factors, - decomposition) + return Emulator{FT}( + machine_learning_tool, + training_pairs, + input_mean, + normalize_inputs, + sqrt_inv_input_cov, + standardize_outputs, + standardize_outputs_factors, + decomposition + ) end """ @@ -158,7 +140,7 @@ end optimize the hyperparameters in the machine learning tool """ -function optimize_hyperparameters!(emulator::Emulator{FT}) where {FT} +function optimize_hyperparameters!(emulator::Emulator{FT}) where {FT<:AbstractFloat} optimize_hyperparameters!(emulator.machine_learning_tool) end @@ -168,8 +150,11 @@ end makes a prediction using the emulator on new inputs (each new inputs given as data columns), default is to predict in the decorrelated space """ -function predict(emulator::Emulator{FT}, new_inputs; transform_to_real=false) where {FT} - +function predict( + emulator::Emulator{FT}, + new_inputs::AbstractMatrix{FT}; + transform_to_real=false +) where {FT<:AbstractFloat} # Check if the size of new_inputs is consistent with the GP model's input # dimension. input_dim, output_dim = size(emulator.training_pairs, 1) @@ -222,7 +207,7 @@ end normalize the input data, with a normalizing function """ -function normalize(emulator::Emulator{FT}, inputs) where {FT} +function normalize(emulator::Emulator{FT}, inputs::AbstractVecOrMat{FT}) where {FT<:AbstractFloat} if emulator.normalize_inputs return normalize(inputs, emulator.input_mean, emulator.sqrt_inv_input_cov) else @@ -235,7 +220,11 @@ end normalize with the empirical Gaussian distribution of points """ -function normalize(inputs, input_mean, sqrt_inv_input_cov) +function normalize( + inputs::AbstractVecOrMat{FT}, + input_mean::AbstractVector{FT}, + sqrt_inv_input_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}} +) where {FT<:AbstractFloat} training_inputs = sqrt_inv_input_cov * (inputs .- input_mean) return training_inputs end @@ -245,33 +234,59 @@ end standardize with a vector of factors (size equal to output dimension) """ -function standardize(outputs, output_cov, factors) - standardized_outputs = outputs ./ factors - standardized_cov = output_cov ./ (factors .* factors') - return standardized_outputs, standardized_cov +function standardize( + outputs::AbstractVecOrMat{FT}, + output_covs::Vector{<:Union{AbstractMatrix{FT}, UniformScaling{FT}}}, + factors::AbstractVector{FT} +) where {FT<:AbstractFloat} + # Case where `output_covs` is a Vector of covariance matrices + outputs = outputs ./ factors + cov_factors = factors .* factors' + for i in 1:length(output_covs) + output_covs[i] = output_covs[i] ./ cov_factors + end + return outputs, output_covs +end + +function standardize( + outputs::AbstractVecOrMat{FT}, + output_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}}, + factors::AbstractVector{FT} +) where {FT<:AbstractFloat} + # Case where `output_cov` is a single covariance matrix + stdized_out, stdized_covs = standardize(outputs, [output_cov], factors) + return stdized_out, stdized_covs[1] end """ - function reverse_standardize(emulator::Emulator, outputs, output_cov) + function reverse_standardize(emulator::Emulator, outputs, output_covs) -reverse a previous standardization with the stored vector of factors (size equal to output dimension) +reverse a previous standardization with the stored vector of factors (size equal to output +dimension). `output_cov` is a Vector of covariance matrices, such as is returned by +[`svd_reverse_transform_mean_cov`](@ref). """ -function reverse_standardize(emulator::Emulator{FT}, outputs, output_cov) where {FT} +function reverse_standardize( + emulator::Emulator{FT}, + outputs::AbstractVecOrMat{FT}, + output_covs::Union{AbstractMatrix{FT}, UniformScaling{FT}, Vector{<:AbstractMatrix{FT}}} +) where {FT<:AbstractFloat} if emulator.standardize_outputs - return standardize(outputs, output_cov, 1.0 ./ emulator.standardize_outputs_factors) + return standardize(outputs, output_covs, 1.0 ./ emulator.standardize_outputs_factors) else - return outputs, output_cov + return outputs, output_covs end end """ -svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, Nothing}) where {FT} + svd_transform(data, obs_noise_cov, truncate_svd) where {FT} Apply a singular value decomposition (SVD) to the data - `data` - GP training data/targets; output_dim × N_samples - `obs_noise_cov` - covariance of observational noise + - `truncate_svd` - Project onto this fraction of the largest principal components. Defaults + to 1.0 (no truncation). Returns the transformed data and the decomposition, which is a matrix factorization of type LinearAlgebra.SVD. @@ -281,102 +296,54 @@ 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; - truncate_svd::FT=1.0) where {FT} - - if obs_noise_cov !== nothing - @assert ndims(obs_noise_cov) == 2 + data::AbstractMatrix{FT}, + obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}, Nothing}; + truncate_svd::FT=1.0 +) where {FT<:AbstractFloat} + if obs_noise_cov === nothing + # no-op case + return data, nothing end - if obs_noise_cov !== nothing + # actually have a matrix to take the SVD of + decomp = svd(obs_noise_cov) + sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomp.S)) + if truncate_svd < 1.0 # Truncate the SVD as a form of regularization - if truncate_svd<1.0 # this variable needs to be provided to this function - # Perform SVD - decomposition = svd(obs_noise_cov) - # Find cutoff - σ = decomposition.S - σ_cumsum = cumsum(σ) / sum(σ); - P_cutoff = truncate_svd; - ind = findall(x->x>P_cutoff, σ_cumsum); k = ind[1] - println("SVD truncated at k:") - println(k) - # Apply truncated SVD - n = size(obs_noise_cov)[1] - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - transformed_data = sqrt_singular_values_inv[1:k,1:k] * decomposition.Vt[1:k,:] * data - transformed_data = transformed_data; - decomposition = Decomposition(decomposition.V[:,1:k], decomposition.Vt[1:k,:], - decomposition.S[1:k], n) + # Find cutoff + S_cumsum = cumsum(decomp.S) / sum(decomp.S) + ind = findall(x -> (x > truncate_svd), S_cumsum) + k = ind[1] + n = size(obs_noise_cov)[1] + println("SVD truncated at k: ", k, "/", n) + # Apply truncated SVD + transformed_data = sqrt_singular_values_inv[1:k, 1:k] * decomp.Vt[1:k, :] * data + decomposition = SVD(decomp.U[:, 1:k], decomp.S[1:k], decomp.Vt[1:k, :]) else - decomposition = svd(obs_noise_cov) - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - transformed_data = sqrt_singular_values_inv * decomposition.Vt * data - decomposition = Decomposition(svd(obs_noise_cov)) - end - else - decomposition = nothing - transformed_data = data + transformed_data = sqrt_singular_values_inv * decomp.Vt * data + decomposition = svd(obs_noise_cov) end - return transformed_data, decomposition end - -""" function svd_transform( - data::Vector{FT}, - obs_noise_cov::Union{Array{FT, 2}, Nothing}; - truncate_svd::FT=1.0) where {FT} - -""" -function svd_transform( - data::Vector{FT}, - obs_noise_cov; - truncate_svd::FT=1.0) where {FT} - - - if obs_noise_cov !== nothing - @assert ndims(obs_noise_cov) == 2 - end - if obs_noise_cov !== nothing - # Truncate the SVD as a form of regularization - if truncate_svd<1.0 # this variable needs to be provided to this function - # Perform SVD - decomposition = svd(obs_noise_cov) - # Find cutoff - σ = decomposition.S - σ_cumsum = cumsum(σ) / sum(σ); - P_cutoff = truncate_svd; - ind = findall(x->x>P_cutoff, σ_cumsum); k = ind[1] - println("SVD truncated at k:") - println(k) - # Apply truncated SVD - n = size(obs_noise_cov)[1] - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - transformed_data = sqrt_singular_values_inv[1:k,1:k] * decomposition.Vt[1:k,:] * data - transformed_data = transformed_data; - decomposition = Decomposition(decomposition.V[:,1:k], decomposition.Vt[1:k,:], - decomposition.S[1:k], n) - else - decomposition = svd(obs_noise_cov) - sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) - transformed_data = sqrt_singular_values_inv * decomposition.Vt * data - decomposition = Decomposition(svd(obs_noise_cov)) - end - else - decomposition = nothing - transformed_data = data - end - - return transformed_data, decomposition + data::AbstractVector{FT}, + obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}, Nothing}; + truncate_svd::FT = 1.0 +) where {FT<:AbstractFloat} + # method for 1D data + transformed_data, decomposition = svd_transform( + reshape(data, :, 1), obs_noise_cov; truncate_svd=truncate_svd + ) + return vec(transformed_data), decomposition end """ -svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::{Array{FT, 2}, decomposition::SVD) where {FT} +svd_reverse_transform_mean_cov(μ, σ2, decomposition::SVD) where {FT} Transform the mean and covariance back to the original (correlated) coordinate system - `μ` - predicted mean; output_dim × N_predicted_points - `σ2` - predicted variance; output_dim × N_predicted_points + - `decomposition` - SVD decomposition of obs_noise_cov. Returns the transformed mean (output_dim × N_predicted_points) and variance. Note that transforming the variance back to the original coordinate system @@ -385,15 +352,16 @@ 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 × output_dim """ -function svd_reverse_transform_mean_cov(μ, σ2, - decomposition::Union{SVD, Decomposition}; - truncate_svd::FT=1.0) where {FT} - @assert ndims(μ) == 2 - @assert ndims(σ2) == 2 +function svd_reverse_transform_mean_cov( + μ::AbstractMatrix{FT}, + σ2::AbstractMatrix{FT}, + decomposition::SVD +) where {FT<:AbstractFloat} + @assert size(μ) == size(σ2) output_dim, N_predicted_points = size(σ2) # We created meanvGP = D_inv * Vt * mean_v so meanv = V * D * meanvGP - sqrt_singular_values= Diagonal(sqrt.(decomposition.S)) + sqrt_singular_values = Diagonal(sqrt.(decomposition.S)) transformed_μ = decomposition.V * sqrt_singular_values * μ transformed_σ2 = [zeros(output_dim, output_dim) for i in 1:N_predicted_points] diff --git a/src/GaussianProcess.jl b/src/GaussianProcess.jl index c2915e6e0..d93c98047 100644 --- a/src/GaussianProcess.jl +++ b/src/GaussianProcess.jl @@ -1,5 +1,7 @@ using GaussianProcesses + +using ..DataStorage using PyCall using ScikitLearn @@ -51,13 +53,13 @@ $(DocStringExtensions.FIELDS) """ struct GaussianProcess{GPPackage} <: MachineLearningTool "the Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs" - models::Vector + models::Vector{Union{<:GaussianProcesses.GPE, <:PyObject, Nothing}} "Kernel object" - kernel + kernel::Union{<:Kernel, <:PyObject, Nothing} "learn the noise with the White Noise kernel explicitly?" - noise_learn + noise_learn::Bool "prediction type (`y` to predict the data, `f` to predict the latent function)" - prediction_type::Union{Nothing, PredictionType} + prediction_type::PredictionType end @@ -75,17 +77,16 @@ end - `prediction_type` - PredictionType object. default predicts data, not latent function (FType()). """ function GaussianProcess( - package; - kernel::Union{K, KPy, Nothing}=nothing, - noise_learn=true, - prediction_type::PredictionType=YType(), - ) where {K<:Kernel, KPy<:PyObject} + package::GPPkg; + kernel::Union{K, KPy, Nothing} = nothing, + noise_learn = true, + prediction_type::PredictionType = YType(), +) where {GPPkg<:GaussianProcessesPackage, K<:Kernel, KPy<:PyObject} # Initialize vector for GP models - models = Any[] + models = Vector{Union{<:GaussianProcesses.GPE, <:PyObject, Nothing}}(undef, 0) return GaussianProcess{typeof(package)}(models, kernel, noise_learn, prediction_type) - end # First we create the GPJL implementation @@ -98,8 +99,8 @@ method to build gaussian process models based on the package """ function build_models!( gp::GaussianProcess{GPJL}, - input_output_pairs) - + input_output_pairs::PairedDataContainer{FT} +) where {FT<:AbstractFloat} # get inputs and outputs input_values = get_inputs(input_output_pairs) output_values = get_outputs(input_output_pairs) @@ -150,7 +151,7 @@ function build_models!( println("kernel in GaussianProcess:") println(kernel_i) data_i = output_values[i,:] - # GPE() arguments: + # GaussianProcesses.GPE() arguments: # input_values: (input_dim × N_samples) # GPdata_i: (N_samples,) @@ -158,7 +159,7 @@ function build_models!( kmean = MeanZero() # Instantiate GP model - m = GPE(input_values, + m = GaussianProcesses.GPE(input_values, output_values[i,:], kmean, kernel_i, @@ -187,52 +188,43 @@ function optimize_hyperparameters!(gp::GaussianProcess{GPJL}) end end -""" - function predict(gp::GaussianProcess{package}, new_inputs::Array{FT, 2}) - -predict means and covariances in decorrelated output space using gaussian process models -""" -predict(gp::GaussianProcess{GPJL}, new_inputs::Array{FT, 2}) where {FT} = predict(gp, new_inputs, gp.prediction_type) - -function predict(gp::GaussianProcess{GPJL}, new_inputs::Array{FT, 2}, ::YType) where {FT} - - +# subroutine with common predict() logic +function _predict( + gp::GaussianProcess, + new_inputs::AbstractMatrix{FT}, + predict_method::Function +) where {FT<:AbstractFloat} M = length(gp.models) - N_samples = size(new_inputs,2) + N_samples = size(new_inputs, 2) # Predicts columns of inputs: input_dim × N_samples μ = zeros(M, N_samples) σ2 = zeros(M, N_samples) for i in 1:M - μ[i,:],σ2[i,:] = predict_y(gp.models[i], new_inputs) + μ[i,:], σ2[i,:] = predict_method(gp.models[i], new_inputs) end - return μ, σ2 end +predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}, ::YType) where {FT<:AbstractFloat} = + _predict(gp, new_inputs, GaussianProcesses.predict_y) +predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}, ::FType) where {FT<:AbstractFloat} = + _predict(gp, new_inputs, GaussianProcesses.predict_f) -function predict(gp::GaussianProcess{GPJL}, new_inputs::Array{FT, 2}, ::FType) where {FT} - - M = length(gp.models) - N_samples = size(new_inputs,2) - # Predicts columns of inputs: input_dim × N_samples - μ = zeros(M, N_samples) - σ2 = zeros(M, N_samples) - for i in 1:M - μ[i,:],σ2[i,:] = predict_f(gp.models[i], new_inputs) - end - - return μ, σ2 -end - +""" + function predict(gp::GaussianProcess{package}, new_inputs) +predict means and covariances in decorrelated output space using gaussian process models +""" +predict(gp::GaussianProcess{GPJL}, new_inputs::AbstractMatrix{FT}) where {FT<:AbstractFloat} = + predict(gp, new_inputs, gp.prediction_type) #now we build the SKLJL implementation function build_models!( gp::GaussianProcess{SKLJL}, - input_output_pairs) - + input_output_pairs::PairedDataContainer{FT} +) where {FT<:AbstractFloat} # get inputs and outputs input_values = permutedims(get_inputs(input_output_pairs), (2,1)) output_values = get_outputs(input_output_pairs) @@ -301,20 +293,11 @@ function optimize_hyperparameters!(gp::GaussianProcess{SKLJL}) println("SKlearn, already trained. continuing...") end - -function predict(gp::GaussianProcess{SKLJL}, new_inputs::Array{FT, 2}) where {FT} - - M = length(gp.models) - N_samples = size(new_inputs,2) - +function _SKJL_predict_function(gp_model::PyObject, new_inputs::AbstractMatrix{FT}) where {FT<:AbstractFloat} # SKJL based on rows not columns; need to transpose inputs - μ = zeros(M, N_samples) - σ = zeros(M, N_samples) - for i in 1:M - μ[i,:],σ[i,:] = gp.models[i].predict(new_inputs', return_std=true) - end - σ2 = σ .* σ - - return μ, σ2 + μ, σ = gp_model.predict(new_inputs', return_std=true) + return μ, (σ .* σ) end +predict(gp::GaussianProcess{SKLJL}, new_inputs::AbstractMatrix{FT}) where {FT<:AbstractFloat} = + _predict(gp, new_inputs, _SKJL_predict_function) diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl index 9f4ef314a..e5a26e26c 100644 --- a/src/MarkovChainMonteCarlo.jl +++ b/src/MarkovChainMonteCarlo.jl @@ -33,23 +33,23 @@ struct MCMC{FT<:AbstractFloat, IT<:Int} "a single sample from the observations. Can e.g. be picked from an Observation struct using get_obs_sample" obs_sample::Vector{FT} "covariance of the observational noise" - obs_noise_cov::Array{FT, 2} + obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}} "array of length N_parameters with the parameters' prior distributions" prior::ParameterDistribution "MCMC step size" - step::Array{FT} + step::FT "Number of MCMC steps that are considered burnin" burnin::IT "the current parameters" - param::Vector{FT} + param::AbstractVector{FT} "Array of accepted MCMC parameter samples. The histogram of these samples gives an approximation of the posterior distribution of the parameters. param_dim x n_samples" - posterior::Array{FT, 2} + posterior::AbstractMatrix{FT} "the (current) value of the logarithm of the posterior (= log_likelihood + log_prior of the current parameters)" - log_posterior::Array{Union{FT, Nothing}} + log_posterior::Union{FT, Nothing} "iteration/step of the MCMC" - iter::Array{IT} + iter::IT "number of accepted proposals" - accept::Array{IT} + accept::IT "MCMC algorithm to use - currently implemented: 'rmw' (random walk Metropolis)" algtype::String "Random number generator object (algorithm + seed) used for sampling and noise, for reproducibility." @@ -70,22 +70,20 @@ where max_iter is the number of MCMC steps to perform (e.g., 100_000) """ function MCMC( - obs_sample::Vector{FT}, - obs_noise_cov::Array{FT, 2}, + obs_sample::AbstractVector{FT}, + obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}}, prior::ParameterDistribution, step::FT, - param_init::Vector{FT}, + param_init::AbstractVector{FT}, max_iter::IT, algtype::String, burnin::IT; - svdflag=true, - standardize=false, - norm_factor::Union{Array{FT, 1}, Nothing}=nothing, - truncate_svd=1.0, - rng::Random.AbstractRNG = Random.GLOBAL_RNG, + svdflag = true, + standardize = false, + norm_factor::Union{AbstractVector{FT}, Nothing} = nothing, + truncate_svd = 1.0, + rng::Random.AbstractRNG = Random.GLOBAL_RNG ) where {FT<:AbstractFloat, IT<:Int} - - param_init_copy = deepcopy(param_init) # Standardize MCMC input? @@ -102,44 +100,45 @@ function MCMC( # We need to transform obs_sample into the correct space if svdflag println("Applying SVD to decorrelating outputs, if not required set svdflag=false") - obs_sample, unused = Emulators.svd_transform(obs_sample, obs_noise_cov; truncate_svd=truncate_svd) + obs_sample, _ = Emulators.svd_transform(obs_sample, obs_noise_cov; truncate_svd=truncate_svd) else println("Assuming independent outputs.") end println(obs_sample) # first row is param_init - posterior = zeros(length(param_init_copy),max_iter + 1) + posterior = zeros(length(param_init_copy), max_iter + 1) posterior[:, 1] = param_init_copy param = param_init_copy - log_posterior = [nothing] - iter = [1] - accept = [0] + log_posterior = nothing + iter = 1 + accept = 0 if algtype != "rwm" error("only random walk metropolis 'rwm' is implemented so far") end - MCMC{FT,IT}(obs_sample, - obs_noise_cov, - prior, - [step], - burnin, - param, - posterior, - log_posterior, - iter, - accept, - algtype, - rng - ) + MCMC{FT,IT}( + obs_sample, + obs_noise_cov, + prior, + step, + burnin, + param, + posterior, + log_posterior, + iter, + accept, + algtype, + rng + ) end -function reset_with_step!(mcmc::MCMC{FT}, step::FT) where {FT} +function reset_with_step!(mcmc::MCMC{FT,IT}, step::FT) where {FT<:AbstractFloat, IT<:Int} # reset to beginning with new stepsize - mcmc.step[1] = step - mcmc.log_posterior[1] = nothing - mcmc.iter[1] = 1 - mcmc.accept[1] = 0 + mcmc.step = step + mcmc.log_posterior = nothing + mcmc.iter = 1 + mcmc.accept = 0 mcmc.posterior[:,2:end] = zeros(size(mcmc.posterior[:,2:end])) mcmc.param[:] = mcmc.posterior[:, 1] end @@ -157,45 +156,54 @@ function get_posterior(mcmc::MCMC) end -function mcmc_sample!(mcmc::MCMC{FT}, g::Vector{FT}, - gcov::Union{Matrix{FT},Diagonal{FT}}) where {FT} +function mcmc_sample!( + mcmc::MCMC{FT,IT}, + g::AbstractVector{FT}, + gcov::Union{AbstractMatrix{FT}, UniformScaling{FT}} +) where {FT<:AbstractFloat, IT<:Int} if mcmc.algtype == "rwm" log_posterior = log_likelihood(mcmc, g, gcov) + log_prior(mcmc) end - if mcmc.log_posterior[1] isa Nothing # do an accept step. - mcmc.log_posterior[1] = log_posterior - log(FT(0.5)) # this makes p_accept = 0.5 + if mcmc.log_posterior === nothing # do an accept step. + mcmc.log_posterior = log_posterior - log(FT(0.5)) # this makes p_accept = 0.5 end # Get new parameters by comparing likelihood_current * prior_current to # likelihood_proposal * prior_proposal - either we accept the proposed # parameter or we stay where we are. - p_accept = exp(log_posterior - mcmc.log_posterior[1]) + p_accept = exp(log_posterior - mcmc.log_posterior) if p_accept > rand(mcmc.rng, Distributions.Uniform(0, 1)) - mcmc.posterior[:,1 + mcmc.iter[1]] = mcmc.param - mcmc.log_posterior[1] = log_posterior - mcmc.accept[1] = mcmc.accept[1] + 1 + mcmc.posterior[:,1 + mcmc.iter] = mcmc.param + mcmc.log_posterior = log_posterior + mcmc.accept = mcmc.accept + 1 else - mcmc.posterior[:,1 + mcmc.iter[1]] = mcmc.posterior[:,mcmc.iter[1]] + mcmc.posterior[:, 1 + mcmc.iter[1]] = mcmc.posterior[:,mcmc.iter] end - mcmc.param[:] = proposal(mcmc)[:] - mcmc.iter[1] = mcmc.iter[1] + 1 + mcmc.param = proposal(mcmc)[:] + mcmc.iter = mcmc.iter + 1 end -function mcmc_sample!(mcmc::MCMC{FT}, g::Vector{FT}, gvar::Vector{FT}) where {FT} +function mcmc_sample!( + mcmc::MCMC{FT,IT}, + g::AbstractVector{FT}, + gvar::AbstractVector{FT} +) where {FT<:AbstractFloat, IT<:Int} return mcmc_sample!(mcmc, g, Diagonal(gvar)) end -function accept_ratio(mcmc::MCMC{FT}) where {FT} - return convert(FT, mcmc.accept[1]) / mcmc.iter[1] +function accept_ratio(mcmc::MCMC{FT, IT}) where {FT<:AbstractFloat, IT<:Int} + return convert(FT, mcmc.accept) / mcmc.iter end -function log_likelihood(mcmc::MCMC{FT}, - g::Vector{FT}, - gcov::Union{Matrix{FT},Diagonal{FT}}) where {FT} - log_rho = FT[0] +function log_likelihood( + mcmc::MCMC{FT,IT}, + g::AbstractVector{FT}, + gcov::Union{AbstractMatrix{FT}, UniformScaling{FT}} +) where {FT<:AbstractFloat, IT<:Int} + log_rho = 0.0 #if gcov == nothing # diff = g - mcmc.obs_sample # log_rho[1] = -FT(0.5) * diff' * (mcmc.obs_noise_cov \ diff) @@ -209,13 +217,13 @@ function log_likelihood(mcmc::MCMC{FT}, log_gpfidelity = -FT(0.5) * sum(log.(eigs)) # Combine got log_rho diff = g - mcmc.obs_sample - log_rho[1] = -FT(0.5) * diff' * (gcov \ diff) + log_gpfidelity + log_rho = -FT(0.5) * diff' * (gcov \ diff) + log_gpfidelity #end - return log_rho[1] + return log_rho end -function log_prior(mcmc::MCMC{FT}) where {FT} +function log_prior(mcmc::MCMC) return get_logpdf(mcmc.prior,mcmc.param) end @@ -228,13 +236,13 @@ function proposal(mcmc::MCMC) prop_dist = MvNormal(zeros(length(mcmc.param)), (mcmc.step[1]^2) * proposal_covariance) end - sample = mcmc.posterior[:,1 + mcmc.iter[1]] .+ rand(mcmc.rng, prop_dist) + sample = mcmc.posterior[:,1 + mcmc.iter] .+ rand(mcmc.rng, prop_dist) return sample end -function find_mcmc_step!(mcmc_test::MCMC{FT}, em::Emulator{FT}; max_iter=2000) where {FT} - step = mcmc_test.step[1] +function find_mcmc_step!(mcmc_test::MCMC{FT,IT}, em::Emulator{FT}; max_iter::IT=2000) where {FT<:AbstractFloat, IT<:Int} + step = mcmc_test.step mcmc_accept = false doubled = false halved = false @@ -291,14 +299,15 @@ function find_mcmc_step!(mcmc_test::MCMC{FT}, em::Emulator{FT}; max_iter=2000) w end - return mcmc_test.step[1] + return mcmc_test.step end -function sample_posterior!(mcmc::MCMC{FT,IT}, - em::Emulator{FT}, - max_iter::IT) where {FT,IT<:Int} - +function sample_posterior!( + mcmc::MCMC{FT, IT}, + em::Emulator{FT}, + max_iter::IT +) where {FT<:AbstractFloat, IT<:Int} for mcmcit in 1:max_iter param = reshape(mcmc.param, :, 1) # test predictions (param is 1 x N_parameters) diff --git a/test/Emulator/runtests.jl b/test/Emulator/runtests.jl index 8610034ae..b50d476ee 100644 --- a/test/Emulator/runtests.jl +++ b/test/Emulator/runtests.jl @@ -33,22 +33,21 @@ using CalibrateEmulateSample.DataContainers # [1.] test SVD (2D y version) test_SVD = svd(Σ) transformed_y, decomposition = Emulators.svd_transform(y, Σ, truncate_svd=1.0) - @test_throws AssertionError Emulators.svd_transform(y, Σ[:,1], truncate_svd=1.0) + @test_throws MethodError Emulators.svd_transform(y, Σ[:,1], truncate_svd=1.0) @test test_SVD.V[:,:] == decomposition.V #(use [:,:] to make it an array) @test test_SVD.Vt == decomposition.Vt @test test_SVD.S == decomposition.S - @test size(test_SVD.S)[1] == decomposition.N @test size(transformed_y) == size(y) # 1D y version transformed_y, decomposition2 = Emulators.svd_transform(y[:, 1], Σ, truncate_svd=1.0) - @test_throws AssertionError Emulators.svd_transform(y[:,1], Σ[:,1], truncate_svd=1.0) + @test_throws MethodError Emulators.svd_transform(y[:,1], Σ[:,1], truncate_svd=1.0) @test size(transformed_y) == size(y[:, 1]) # Reverse SVD - y_new, y_cov_new = Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y,(d,1)),ones(d,1),decomposition2, truncate_svd=1.0) - @test_throws AssertionError Emulators.svd_reverse_transform_mean_cov(transformed_y,ones(d,1),decomposition2, truncate_svd=1.0) - @test_throws AssertionError Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y,(d,1)),ones(d),decomposition2, truncate_svd=1.0) + y_new, y_cov_new = Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y, (d,1)), ones(d,1), decomposition2) + @test_throws MethodError Emulators.svd_reverse_transform_mean_cov(transformed_y, ones(d,1), decomposition2) + @test_throws MethodError Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y,(d,1)), ones(d), decomposition2) @test size(y_new)[1] == size(y[:, 1])[1] @test y_new ≈ y[:,1] @test y_cov_new[1] ≈ Σ @@ -60,7 +59,7 @@ using CalibrateEmulateSample.DataContainers @test size(transformed_y)[1] == trunc_size # [2.] test Normalization - input_mean = reshape(mean(get_inputs(iopairs), dims=2), :, 1) #column vector + input_mean = vec(mean(get_inputs(iopairs), dims=2)) #column vector sqrt_inv_input_cov = sqrt(inv(Symmetric(cov(get_inputs(iopairs), dims=2)))) norm_inputs = Emulators.normalize(get_inputs(iopairs),input_mean,sqrt_inv_input_cov) @@ -107,7 +106,6 @@ using CalibrateEmulateSample.DataContainers @test test_decomp.V == decomposition.V #(use [:,:] to make it an array) @test test_decomp.Vt == decomposition.Vt @test test_decomp.S == decomposition.S - @test test_decomp.N == decomposition.N gp = GaussianProcess(GPJL()) emulator2 = Emulator( From 8c0550d692e04ce0c226f6ab3862cf3f9362d4c6 Mon Sep 17 00:00:00 2001 From: Thomas Jackson Date: Mon, 25 Apr 2022 23:28:34 -0400 Subject: [PATCH 2/4] Fix DataContainers import --- src/GaussianProcess.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/GaussianProcess.jl b/src/GaussianProcess.jl index d93c98047..afb4b2880 100644 --- a/src/GaussianProcess.jl +++ b/src/GaussianProcess.jl @@ -1,7 +1,7 @@ using GaussianProcesses -using ..DataStorage +using EnsembleKalmanProcesses.DataContainers using PyCall using ScikitLearn From f221db27159de8fe991847d2c2f317f637fee1a8 Mon Sep 17 00:00:00 2001 From: Thomas Jackson Date: Mon, 25 Apr 2022 23:36:55 -0400 Subject: [PATCH 3/4] Make MCMC struct mutable (to be reverted in PR #130) --- src/MarkovChainMonteCarlo.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl index e5a26e26c..5af853ab1 100644 --- a/src/MarkovChainMonteCarlo.jl +++ b/src/MarkovChainMonteCarlo.jl @@ -29,7 +29,7 @@ Structure to organize MCMC parameters and data # Fields $(DocStringExtensions.FIELDS) """ -struct MCMC{FT<:AbstractFloat, IT<:Int} +mutable struct MCMC{FT<:AbstractFloat, IT<:Int} "a single sample from the observations. Can e.g. be picked from an Observation struct using get_obs_sample" obs_sample::Vector{FT} "covariance of the observational noise" From cbaf482a660e9b57add55f4918fef32f1901fee6 Mon Sep 17 00:00:00 2001 From: Thomas Jackson Date: Tue, 26 Apr 2022 01:39:34 -0400 Subject: [PATCH 4/4] Methods & tests for Diagonal, UniformScaling cov types --- src/Emulator.jl | 38 +++++++- test/Emulator/runtests.jl | 190 +++++++++++++++++++++----------------- 2 files changed, 139 insertions(+), 89 deletions(-) diff --git a/src/Emulator.jl b/src/Emulator.jl index d3fca75ea..8b6f12873 100644 --- a/src/Emulator.jl +++ b/src/Emulator.jl @@ -73,6 +73,11 @@ function Emulator( # For Consistency checks input_dim, output_dim = size(input_output_pairs, 1) + if isa(obs_noise_cov, UniformScaling) + # Cast UniformScaling to Diagonal, since UniformScaling incompatible with + # SVD and standardize() + obs_noise_cov = Diagonal(obs_noise_cov, output_dim) + end 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)) @@ -240,10 +245,25 @@ function standardize( factors::AbstractVector{FT} ) where {FT<:AbstractFloat} # Case where `output_covs` is a Vector of covariance matrices + n = length(factors) outputs = outputs ./ factors cov_factors = factors .* factors' for i in 1:length(output_covs) - output_covs[i] = output_covs[i] ./ cov_factors + if isa(output_covs[i], Matrix) + output_covs[i] = output_covs[i] ./ cov_factors + elseif isa(output_covs[i], UniformScaling) + # assert all elements of factors are same, otherwise can't cast back to + # UniformScaling + @assert all(y -> y == first(factors), factors) + output_covs[i] = output_covs[i] / cov_factors[1,1] + else + # Diagonal, TriDiagonal etc. case + # need to do ./ as dense matrix and cast back to original type + # https://discourse.julialang.org/t/get-generic-constructor-of-parametric-type/57189/8 + T = typeof(output_covs[i]) + cast_T = getfield(parentmodule(T), nameof(T)) + output_covs[i] = cast_T(Matrix(output_covs[i]) ./ cov_factors) + end end return outputs, output_covs end @@ -268,7 +288,7 @@ dimension). `output_cov` is a Vector of covariance matrices, such as is returned function reverse_standardize( emulator::Emulator{FT}, outputs::AbstractVecOrMat{FT}, - output_covs::Union{AbstractMatrix{FT}, UniformScaling{FT}, Vector{<:AbstractMatrix{FT}}} + output_covs::Union{AbstractMatrix{FT}, Vector{<:AbstractMatrix{FT}}} ) where {FT<:AbstractFloat} if emulator.standardize_outputs return standardize(outputs, output_covs, 1.0 ./ emulator.standardize_outputs_factors) @@ -297,7 +317,7 @@ in S are sorted in descending order. """ function svd_transform( data::AbstractMatrix{FT}, - obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}, Nothing}; + obs_noise_cov::Union{AbstractMatrix{FT}, Nothing}; truncate_svd::FT=1.0 ) where {FT<:AbstractFloat} if obs_noise_cov === nothing @@ -313,7 +333,7 @@ function svd_transform( S_cumsum = cumsum(decomp.S) / sum(decomp.S) ind = findall(x -> (x > truncate_svd), S_cumsum) k = ind[1] - n = size(obs_noise_cov)[1] + n = size(data)[1] println("SVD truncated at k: ", k, "/", n) # Apply truncated SVD transformed_data = sqrt_singular_values_inv[1:k, 1:k] * decomp.Vt[1:k, :] * data @@ -327,7 +347,7 @@ end function svd_transform( data::AbstractVector{FT}, - obs_noise_cov::Union{AbstractMatrix{FT}, UniformScaling{FT}, Nothing}; + obs_noise_cov::Union{AbstractMatrix{FT}, Nothing}; truncate_svd::FT = 1.0 ) where {FT<:AbstractFloat} # method for 1D data @@ -337,6 +357,14 @@ function svd_transform( return vec(transformed_data), decomposition end +function svd_transform( + data::AbstractVecOrMat{FT}, obs_noise_cov::UniformScaling{FT}; + truncate_svd::FT = 1.0 +) where {FT<:AbstractFloat} + # method for UniformScaling + return svd_transform(data, Diagonal(obs_noise_cov, size(data, 1)); truncate_svd=truncate_svd) +end + """ svd_reverse_transform_mean_cov(μ, σ2, decomposition::SVD) where {FT} diff --git a/test/Emulator/runtests.jl b/test/Emulator/runtests.jl index b50d476ee..ef35e5dee 100644 --- a/test/Emulator/runtests.jl +++ b/test/Emulator/runtests.jl @@ -8,80 +8,22 @@ using LinearAlgebra using CalibrateEmulateSample.Emulators using CalibrateEmulateSample.DataContainers -@testset "Emulators" begin - - - #build some quick data + noise - m=50 - d=6 - x = rand(3, m) #R^3 - y = rand(d, m) #R^5 - - # "noise" - μ = zeros(d) - Σ = rand(d,d) - Σ = Σ'*Σ - noise_samples = rand(MvNormal(μ, Σ), m) - y += noise_samples - - iopairs = PairedDataContainer(x,y,data_are_columns=true) - @test get_inputs(iopairs) == x - @test get_outputs(iopairs) == y - - - - # [1.] test SVD (2D y version) - test_SVD = svd(Σ) - transformed_y, decomposition = Emulators.svd_transform(y, Σ, truncate_svd=1.0) - @test_throws MethodError Emulators.svd_transform(y, Σ[:,1], truncate_svd=1.0) - @test test_SVD.V[:,:] == decomposition.V #(use [:,:] to make it an array) - @test test_SVD.Vt == decomposition.Vt - @test test_SVD.S == decomposition.S - @test size(transformed_y) == size(y) - - # 1D y version - transformed_y, decomposition2 = Emulators.svd_transform(y[:, 1], Σ, truncate_svd=1.0) - @test_throws MethodError Emulators.svd_transform(y[:,1], Σ[:,1], truncate_svd=1.0) - @test size(transformed_y) == size(y[:, 1]) - - # Reverse SVD - y_new, y_cov_new = Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y, (d,1)), ones(d,1), decomposition2) - @test_throws MethodError Emulators.svd_reverse_transform_mean_cov(transformed_y, ones(d,1), decomposition2) - @test_throws MethodError Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y,(d,1)), ones(d), decomposition2) - @test size(y_new)[1] == size(y[:, 1])[1] - @test y_new ≈ y[:,1] - @test y_cov_new[1] ≈ Σ - - # Truncation - transformed_y, trunc_decomposition = Emulators.svd_transform(y[:, 1], Σ, truncate_svd=0.95) - trunc_size = size(trunc_decomposition.S)[1] - @test test_SVD.S[1:trunc_size] == trunc_decomposition.S - @test size(transformed_y)[1] == trunc_size - +#build an unknown type +struct MLTester <: Emulators.MachineLearningTool end + +function constructor_tests( + iopairs::PairedDataContainer, + Σ::Union{AbstractMatrix{FT}, UniformScaling{FT}, Nothing}, + norm_factors, decomposition +) where {FT<:AbstractFloat} # [2.] test Normalization input_mean = vec(mean(get_inputs(iopairs), dims=2)) #column vector sqrt_inv_input_cov = sqrt(inv(Symmetric(cov(get_inputs(iopairs), dims=2)))) - - norm_inputs = Emulators.normalize(get_inputs(iopairs),input_mean,sqrt_inv_input_cov) + norm_inputs = Emulators.normalize(get_inputs(iopairs), input_mean, sqrt_inv_input_cov) @test norm_inputs == sqrt_inv_input_cov * (get_inputs(iopairs) .- input_mean) - # [3.] test Standardization - norm_factors = 10.0 - norm_factors = fill(norm_factors, size(y[:,1])) # must be size of output dim - - s_y, s_y_cov = Emulators.standardize(get_outputs(iopairs),Σ,norm_factors) - @test s_y == get_outputs(iopairs)./norm_factors - @test s_y_cov == Σ ./ (norm_factors.*norm_factors') - - - # [4.] test emulator preserves the structures - - #build an unknown type - struct MLTester <: Emulators.MachineLearningTool end - mlt = MLTester() - @test_throws ErrorException emulator = Emulator( mlt, iopairs, @@ -92,7 +34,6 @@ using CalibrateEmulateSample.DataContainers #build a known type, with defaults gp = GaussianProcess(GPJL()) - emulator = Emulator( gp, iopairs, @@ -116,10 +57,9 @@ using CalibrateEmulateSample.DataContainers standardize_outputs=false, truncate_svd=1.0) train_inputs = get_inputs(emulator2.training_pairs) - @test norm_inputs == train_inputs - - train_inputs2 = Emulators.normalize(emulator2,get_inputs(iopairs)) - @test norm_inputs == train_inputs + @test train_inputs == norm_inputs + train_inputs2 = Emulators.normalize(emulator2, get_inputs(iopairs)) + @test train_inputs2 == norm_inputs # reverse standardise gp = GaussianProcess(GPJL()) @@ -133,24 +73,106 @@ using CalibrateEmulateSample.DataContainers truncate_svd=1.0) #standardized and decorrelated (sd) data + s_y, s_y_cov = Emulators.standardize(get_outputs(iopairs), Σ, norm_factors) sd_train_outputs = get_outputs(emulator3.training_pairs) sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(emulator3.decomposition.S)) decorrelated_s_y = sqrt_singular_values_inv * emulator3.decomposition.Vt * s_y @test decorrelated_s_y == sd_train_outputs +end +@testset "Emulators" begin + #build some quick data + noise + m=50 + d=6 + x = rand(3, m) #R^3 + y = rand(d, m) #R^5 + + # "noise" + μ = zeros(d) + Σ = rand(d,d) + Σ = Σ'*Σ + noise_samples = rand(MvNormal(μ, Σ), m) + y += noise_samples + + iopairs = PairedDataContainer(x,y,data_are_columns=true) + @test get_inputs(iopairs) == x + @test get_outputs(iopairs) == y + + # [1.] test SVD (2D y version) + test_SVD = svd(Σ) + transformed_y, decomposition = Emulators.svd_transform(y, Σ, truncate_svd=1.0) + @test_throws MethodError Emulators.svd_transform(y, Σ[:,1], truncate_svd=1.0) + @test test_SVD.V[:,:] == decomposition.V #(use [:,:] to make it an array) + @test test_SVD.Vt == decomposition.Vt + @test test_SVD.S == decomposition.S + @test size(transformed_y) == size(y) - # truncation - gp = GaussianProcess(GPJL()) - emulator4 = Emulator( - gp, - iopairs, - obs_noise_cov=Σ, - normalize_inputs=false, - standardize_outputs=false, - truncate_svd=0.9) - trunc_size = size(emulator4.decomposition.S)[1] - @test test_SVD.S[1:trunc_size] == emulator4.decomposition.S - + # 1D y version + transformed_y, decomposition2 = Emulators.svd_transform(y[:, 1], Σ, truncate_svd=1.0) + @test_throws MethodError Emulators.svd_transform(y[:,1], Σ[:,1], truncate_svd=1.0) + @test size(transformed_y) == size(y[:, 1]) + # Reverse SVD + y_new, y_cov_new = Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y, (d,1)), ones(d,1), decomposition2) + @test_throws MethodError Emulators.svd_reverse_transform_mean_cov(transformed_y, ones(d,1), decomposition2) + @test_throws MethodError Emulators.svd_reverse_transform_mean_cov(reshape(transformed_y,(d,1)), ones(d), decomposition2) + @test size(y_new)[1] == size(y[:, 1])[1] + @test y_new ≈ y[:,1] + @test y_cov_new[1] ≈ Σ + + # Truncation + transformed_y, trunc_decomposition = Emulators.svd_transform(y[:, 1], Σ, truncate_svd=0.95) + trunc_size = size(trunc_decomposition.S)[1] + @test test_SVD.S[1:trunc_size] == trunc_decomposition.S + @test size(transformed_y)[1] == trunc_size + + # [3.] test Standardization + norm_factors = 10.0 + norm_factors = fill(norm_factors, size(y[:,1])) # must be size of output dim + s_y, s_y_cov = Emulators.standardize(get_outputs(iopairs),Σ,norm_factors) + @test s_y == get_outputs(iopairs)./norm_factors + @test s_y_cov == Σ ./ (norm_factors.*norm_factors') + + @testset "Emulators - dense Array Σ" begin + constructor_tests(iopairs, Σ, norm_factors, decomposition) + + # truncation - only test here, where singular value of Σ differ + gp = GaussianProcess(GPJL()) + emulator4 = Emulator( + gp, + iopairs, + obs_noise_cov=Σ, + normalize_inputs=false, + standardize_outputs=false, + truncate_svd=0.9) + trunc_size = size(emulator4.decomposition.S)[1] + @test test_SVD.S[1:trunc_size] == emulator4.decomposition.S + end + + @testset "Emulators - Diagonal Σ" begin + x = rand(3, m) #R^3 + y = rand(d, m) #R^5 + # "noise" + Σ = Diagonal(rand(d) .+ 1.0) + noise_samples = rand(MvNormal(zeros(d), Σ), m) + y += noise_samples + iopairs = PairedDataContainer(x, y, data_are_columns=true) + _, decomposition = Emulators.svd_transform(y, Σ, truncate_svd=1.0) + + constructor_tests(iopairs, Σ, norm_factors, decomposition) + end + + @testset "Emulators - UniformScaling Σ" begin + x = rand(3, m) #R^3 + y = rand(d, m) #R^5 + # "noise" + Σ = UniformScaling(1.3) + noise_samples = rand(MvNormal(zeros(d), Σ), m) + y += noise_samples + iopairs = PairedDataContainer(x, y, data_are_columns=true) + _, decomposition = Emulators.svd_transform(y, Σ, truncate_svd=1.0) + + constructor_tests(iopairs, Σ, norm_factors, decomposition) + end end