Skip to content

Commit

Permalink
format
Browse files Browse the repository at this point in the history
  • Loading branch information
odunbar committed Oct 20, 2023
1 parent 91bd7a0 commit fc23a09
Show file tree
Hide file tree
Showing 4 changed files with 57 additions and 60 deletions.
98 changes: 47 additions & 51 deletions examples/Emulator/Ishigami/emulate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,34 +28,30 @@ function main()
n_repeats = 20 # repeat exp with same data.

# To create the sampling
n_data_gen = 2000
n_data_gen = 2000

data = SobolData(
params = OrderedDict(
:x1 => Uniform(-π, π),
:x2 => Uniform(-π, π),
:x3 => Uniform(-π, π),
),
N = n_data_gen
params = OrderedDict(:x1 => Uniform(-π, π), :x2 => Uniform(-π, π), :x3 => Uniform(-π, π)),
N = n_data_gen,
)

# To perform global analysis,
# one must generate samples using Sobol sequence (i.e. creates more than N points)
samples = GSA.sample(data)
n_data = size(samples,1) # [n_samples x 3]
n_data = size(samples, 1) # [n_samples x 3]
# run model (example)
y = GSA.ishigami(samples)
# perform Sobol Analysis
result = analyze(data, y);
f1 = Figure(resolution = (1.618*900, 300), markersize=4)
result = analyze(data, y)

f1 = Figure(resolution = (1.618 * 900, 300), markersize = 4)
axx = Axis(f1[1, 1], xlabel = "x1", ylabel = "f")
axy = Axis(f1[1, 2], xlabel = "x2", ylabel = "f")
axz = Axis(f1[1, 3], xlabel = "x3", ylabel = "f")
scatter!(axx, samples[:,1], y[:], color = :orange)
scatter!(axy, samples[:,2], y[:], color = :orange)
scatter!(axz, samples[:,3], y[:], color = :orange)

scatter!(axx, samples[:, 1], y[:], color = :orange)
scatter!(axy, samples[:, 2], y[:], color = :orange)
scatter!(axz, samples[:, 3], y[:], color = :orange)

save("ishigami_slices_truth.png", f1, px_per_unit = 3)
save("ishigami_slices_truth.pdf", f1, px_per_unit = 3)
Expand All @@ -73,7 +69,7 @@ function main()
output[i] = y[ind[i]] + noise[i]
end
iopairs = PairedDataContainer(input, output)

cases = ["Prior", "GP", "RF-scalar"]
case = cases[3]
decorrelate = true
Expand All @@ -94,7 +90,7 @@ function main()
y_preds = []
result_preds = []

for rep_idx = 1:n_repeats
for rep_idx in 1:n_repeats

# Build ML tools
if case == "GP"
Expand All @@ -106,10 +102,10 @@ function main()
prediction_type = pred_type,
noise_learn = false,
)

elseif case ["RF-scalar", "Prior"]
kernel_structure = SeparableKernel(LowRankFactor(3, nugget),OneDimFactor())

kernel_structure = SeparableKernel(LowRankFactor(3, nugget), OneDimFactor())
n_features = 500
mlt = ScalarRandomFeatureInterface(
n_features,
Expand All @@ -119,38 +115,38 @@ function main()
optimizer_options = overrides,
)
end

# Emulate
emulator = Emulator(mlt, iopairs; obs_noise_cov = Γ*I, decorrelate = decorrelate)
emulator = Emulator(mlt, iopairs; obs_noise_cov = Γ * I, decorrelate = decorrelate)
optimize_hyperparameters!(emulator)

# predict on all Sobol points with emulator (example)
y_pred, y_var = predict(emulator, samples', transform_to_real = true)

# obtain emulated Sobol indices
result_pred = analyze(data, y_pred');
result_pred = analyze(data, y_pred')
push!(y_preds, y_pred)
push!(result_preds, result_pred)

end

# analytic sobol indices
a = 7
b = 0.1
V = a^2/8 + b * π^4/5 + b^2 * π^8/18 + 1/2
V1 = 0.5 * (1 + b*π^4/5)^2
V2 = a^2/8
V = a^2 / 8 + b * π^4 / 5 + b^2 * π^8 / 18 + 1 / 2
V1 = 0.5 * (1 + b * π^4 / 5)^2
V2 = a^2 / 8
V3 = 0
VT1 = 0.5*(1+b*π^4/5)^2 + 8*b^2*π^8/225
VT2 = a^2/8
VT3 = 8*b^2*π^8/225
VT1 = 0.5 * (1 + b * π^4 / 5)^2 + 8 * b^2 * π^8 / 225
VT2 = a^2 / 8
VT3 = 8 * b^2 * π^8 / 225



println(" ")
println("True Sobol Indices")
println("******************")
println(" firstorder: ", [V1/V, V2/V, V3/V])
println(" totalorder: ", [VT1/V, VT2/V, VT3/V])
println(" firstorder: ", [V1 / V, V2 / V, V3 / V])
println(" totalorder: ", [VT1 / V, VT2 / V, VT3 / V])
println(" ")
println("Sampled truth Sobol Indices (# points $n_data)")
println("***************************")
Expand All @@ -161,38 +157,38 @@ function main()
println("Sampled Emulated Sobol Indices (# obs $n_train_pts, noise var )")
println("***************************************************************")
if n_repeats == 1
println(" firstorder: ", result_preds[1][:firstorder])
println(" firstorder: ", result_preds[1][:firstorder])
println(" totalorder: ", result_preds[1][:totalorder])
else
else
firstorder_mean = mean([rp[:firstorder] for rp in result_preds])
firstorder_std = std([rp[:firstorder] for rp in result_preds])
totalorder_mean = mean([rp[:totalorder] for rp in result_preds])
totalorder_std = std([rp[:totalorder] for rp in result_preds])

println("(mean) firstorder: ", firstorder_mean)
println("(std) firstorder: ", firstorder_std)
println("(mean) totalorder: ", totalorder_mean)
println("(std) totalorder: ", totalorder_std)
end
end



# plots
f2 = Figure(resolution = (1.618*900, 300), markersize=4)

f2 = Figure(resolution = (1.618 * 900, 300), markersize = 4)
axx_em = Axis(f2[1, 1], xlabel = "x1", ylabel = "f")
axy_em = Axis(f2[1, 2], xlabel = "x2", ylabel = "f")
axz_em = Axis(f2[1, 3], xlabel = "x3", ylabel = "f")
scatter!(axx_em, samples[:,1], y_preds[1][:], color = :blue)
scatter!(axy_em, samples[:,2], y_preds[1][:], color = :blue)
scatter!(axz_em, samples[:,3], y_preds[1][:], color = :blue)
scatter!(axx_em, samples[ind,1], y[ind] + noise, color = :red, markersize=8)
scatter!(axy_em, samples[ind,2], y[ind] + noise, color = :red, markersize=8)
scatter!(axz_em, samples[ind,3], y[ind] + noise, color = :red, markersize=8)
scatter!(axx_em, samples[:, 1], y_preds[1][:], color = :blue)
scatter!(axy_em, samples[:, 2], y_preds[1][:], color = :blue)
scatter!(axz_em, samples[:, 3], y_preds[1][:], color = :blue)
scatter!(axx_em, samples[ind, 1], y[ind] + noise, color = :red, markersize = 8)
scatter!(axy_em, samples[ind, 2], y[ind] + noise, color = :red, markersize = 8)
scatter!(axz_em, samples[ind, 3], y[ind] + noise, color = :red, markersize = 8)

save("ishigami_slices_$(case).png", f2, px_per_unit = 3)
save("ishigami_slices_$(case).pdf", f2, px_per_unit = 3)


end


Expand Down
15 changes: 8 additions & 7 deletions examples/Emulator/L63/emulate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,8 @@ function main()
iopairs = PairedDataContainer(input, output)

# Emulate
cases = ["GP", "RF-scalar", "RF-scalar-diagin", "RF-vector-svd-nonsep","RF-vector-nosvd-nonsep", "RF-vector-nosvd-sep"]
cases =
["GP", "RF-scalar", "RF-scalar-diagin", "RF-vector-svd-nonsep", "RF-vector-nosvd-nonsep", "RF-vector-nosvd-sep"]

case = cases[1]
decorrelate = true
Expand All @@ -74,8 +75,8 @@ function main()
"cov_sample_multiplier" => 0.5,
"n_features_opt" => 200,
"n_iteration" => 20,
# "n_ensemble" => 20,
# "localization" => EKP.Localizers.SEC(1.0,0.1),
# "n_ensemble" => 20,
# "localization" => EKP.Localizers.SEC(1.0,0.1),
"accelerator" => NesterovAccelerator(),
)

Expand Down Expand Up @@ -117,7 +118,7 @@ function main()
elseif case ["RF-vector-nosvd-nonsep"]
kernel_structure = NonseparableKernel(LowRankFactor(3, nugget))
n_features = 500
decorrelate=false # don't do SVD
decorrelate = false # don't do SVD
mlt = VectorRandomFeatureInterface(
n_features,
3,
Expand All @@ -127,9 +128,9 @@ function main()
optimizer_options = overrides,
)
elseif case ["RF-vector-nosvd-sep"]
kernel_structure = SeparableKernel(LowRankFactor(3, nugget),LowRankFactor(3,nugget) )
kernel_structure = SeparableKernel(LowRankFactor(3, nugget), LowRankFactor(3, nugget))
n_features = 500
decorrelate=false # don't do SVD
decorrelate = false # don't do SVD
mlt = VectorRandomFeatureInterface(
n_features,
3,
Expand All @@ -139,7 +140,7 @@ function main()
optimizer_options = overrides,
)
end

# Emulate
emulator = Emulator(mlt, iopairs; obs_noise_cov = Γy, decorrelate = decorrelate)
optimize_hyperparameters!(emulator)
Expand Down
2 changes: 1 addition & 1 deletion src/ScalarRandomFeature.jl
Original file line number Diff line number Diff line change
Expand Up @@ -312,7 +312,7 @@ function build_models!(
decomp_type = get_feature_decomposition(srfi)
optimizer_options = get_optimizer_options(srfi)


# Optimize features with EKP for each output dim
# [1.] Split data into test/train 80/20
train_fraction = optimizer_options["train_fraction"]
Expand Down
2 changes: 1 addition & 1 deletion src/VectorRandomFeature.jl
Original file line number Diff line number Diff line change
Expand Up @@ -505,7 +505,7 @@ function build_models!(
localization = optimizer_options["localization"]
accelerator = optimizer_options["accelerator"]

if !isposdef(Γ)
if !isposdef(Γ)
Γ = posdef_correct(Γ)
end

Expand Down

0 comments on commit fc23a09

Please sign in to comment.