Skip to content

Commit

Permalink
Methods & tests for Diagonal, UniformScaling cov types
Browse files Browse the repository at this point in the history
  • Loading branch information
@tsj5
tsj5 committed Apr 26, 2022
1 parent f221db2 commit cbaf482
Show file tree
Hide file tree
Showing 2 changed files with 139 additions and 89 deletions.
38 changes: 33 additions & 5 deletions src/Emulator.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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}
Expand Down
190 changes: 106 additions & 84 deletions test/Emulator/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -92,7 +34,6 @@ using CalibrateEmulateSample.DataContainers

#build a known type, with defaults
gp = GaussianProcess(GPJL())

emulator = Emulator(
gp,
iopairs,
Expand All @@ -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())
Expand All @@ -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

0 comments on commit cbaf482

Please sign in to comment.