Skip to content

Commit a318edb

Browse files
Replacement of Emulator processing with new encoding framework (#367)
* initial gutting of Emulator transformation and replace with new encoders * emulator bug-fixes * add new interface * remove old normalization/standardizations code * works when no scheduler provided * trunc -> trunc_val * update and cleanup example and plots * gp case * cycle over cases and encoders * n=200 * smaller noise and new schedule train * more flexible regularizer (now it is not always going to be "I" * rm unnecessary plots from compare_regression * new message * converted Ishigami * user_... to ... * emulate * copy -> deepcopy * G-function update * utils * L63 * L63 * sensible default for EKI-produced data * new encode/decode works with sampling too * sinusoid * new encoder methods and exports * bugfix truncation of variance * use new method to encode * fix regularization for non-I matrices * update Lorenz and Lorenz spatial dep examples * updated darcy * updated cloudy example * remove standardization * updated EDMF-data * updated GCM * typo * typo * added more API and == * add ! and mor eunit test coverage for emulators * remove large comment * more cases * GP test updated * use def. enc. for GP tests * all return types * output cases * updated RF tests * update interfcace to abstractGP * updated run tests for MCMC * format * missed by foramtter * docs * docs * docstrings and comments addressed * format * add docs page * add docs into index * Reduce number of methods; give data processors separate files * Pass both structure matrices to PDCPs * Format code * Fix bugs * Adapt tests to new functions * replace logic for input obs * format and emulate docs page * n_iter scalar * n_iter -> scalar * remove default structure matrix, extend framework of encoding a nothing, add tests * behaviour when passing nothing into ML tools * format * when not learning noise, regularization is scaled with output matrix * remove @info * typo * missing \lambda in uniform scaling case * typo * format * missed format file --------- Co-authored-by: Arne Bouillon <arne.bouillon@kuleuven.be> Co-authored-by: ArneBouillon <45404227+ArneBouillon@users.noreply.github.com>
1 parent 212c772 commit a318edb

47 files changed

Lines changed: 1808 additions & 2407 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

docs/make.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,7 @@ pages = [
4242
"Calibrate" => "calibrate.md",
4343
"Emulate" => emulate,
4444
"Sample" => "sample.md",
45+
"Data Processing and Dimension reduction" => "data_processing.md",
4546
"Glossary" => "glossary.md",
4647
"API" => api,
4748
]

docs/src/API/Emulators.md

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -9,9 +9,8 @@ Emulator
99
optimize_hyperparameters!(::Emulator)
1010
Emulator(::MachineLearningTool, ::PairedDataContainer{FT}) where {FT <: AbstractFloat}
1111
predict
12-
normalize
13-
standardize
14-
reverse_standardize
15-
svd_transform
16-
svd_reverse_transform_mean_cov
12+
encode_data
13+
decode_data
14+
encode_structure_matrix
15+
decode_structure_matrix
1716
```

docs/src/API/GaussianProcess.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ GaussianProcess(
1515
::FT,
1616
::PredictionType,
1717
) where {GPPkg <: GaussianProcessesPackage, K <: GaussianProcesses.Kernel, KPy <: PyObject, AK <:AbstractGPs.Kernel, FT <: AbstractFloat}
18-
build_models!(::GaussianProcess{GPJL}, ::PairedDataContainer{FT}) where {FT <: AbstractFloat}
18+
build_models!(::GaussianProcess{GPJL}, ::PairedDataContainer{FT}, input_structure_matrix, output_structure_matrix) where {FT <: AbstractFloat}
1919
optimize_hyperparameters!(::GaussianProcess{GPJL})
2020
predict(::GaussianProcess{GPJL}, ::AbstractMatrix{FT}) where {FT <: AbstractFloat}
2121
```

docs/src/API/MarkovChainMonteCarlo.md

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -50,9 +50,3 @@ EmulatorPosteriorModel
5050
MCMCState
5151
accept_ratio
5252
```
53-
54-
## Internals - Other
55-
56-
```@docs
57-
to_decorrelated
58-
```

docs/src/API/RandomFeatures.md

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ build_default_prior
2222
```@docs
2323
ScalarRandomFeatureInterface
2424
ScalarRandomFeatureInterface(::Int,::Int)
25-
build_models!(::ScalarRandomFeatureInterface, ::PairedDataContainer{FT}) where {FT <: AbstractFloat}
25+
build_models!(::ScalarRandomFeatureInterface, ::PairedDataContainer{FT}, input_structure_matrix, output_structure_matrix) where {FT <: AbstractFloat}
2626
predict(::ScalarRandomFeatureInterface, ::M) where {M <: AbstractMatrix}
2727
```
2828

@@ -31,7 +31,7 @@ predict(::ScalarRandomFeatureInterface, ::M) where {M <: AbstractMatrix}
3131
```@docs
3232
VectorRandomFeatureInterface
3333
VectorRandomFeatureInterface(::Int, ::Int, ::Int)
34-
build_models!(::VectorRandomFeatureInterface, ::PairedDataContainer{FT}) where {FT <: AbstractFloat}
34+
build_models!(::VectorRandomFeatureInterface, ::PairedDataContainer{FT}, input_structure_matrix, output_structure_matrix) where {FT <: AbstractFloat}
3535
predict(::VectorRandomFeatureInterface, ::M) where {M <: AbstractMatrix}
3636
```
3737

docs/src/data_processing.md

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
# [Data Processing and Dimension Reduction](@id data-proc)
2+
3+
## Overview
4+
When working with high-dimensional problems with modest training data pairs, the bottleneck of CES procedure is the training of a competent emulator. It is often necessary to process and dimensionally-reduce the data to ease the learning task of the emulator. We provide a flexible (and extensible!) framework to create encoders and decoders for this purpose. The framework works as follows
5+
- An `encoder_schedule` defines the type of processing to be applied to input and/or output spaces
6+
- The `encoder_schedule` is passed into the `Emulator` where it is initialized and stored.
7+
- The encoder will be used automatically to encode training data and predictions within the Emulate and Sample routines.
8+
9+
An external API is also available using the `encode_data`, and `encode_structure_matrix` methods if needed.
10+
11+
## Define an encoder schedule
12+
13+
The user may provide an encoder schedule to transform data in a useful way. For example, mapping all output data each dimension to be bounded in [0,1]
14+
```julia
15+
simple_schedule = (minmax_scale(), "out")
16+
```
17+
The unit of the encoder contains a `DataProcessor`, the type of processing, and an string, whether to apply to input data,`"in"`, or output data `"out"`, or both `"in_and_out"`.
18+
19+
The encoder schedule can be a vector of several units that apply multiple `DataProcessors` in order:
20+
```julia
21+
complex_schedule = [
22+
(decorrelate_sample_cov(), "in"),
23+
(quartile_scale(), "in"),
24+
(decorrelate_structure_mat(retain_var=0.95), "out"),
25+
(canonical_correlation(), "in_and_out"),
26+
]
27+
```
28+
In this (rather unrealistic) chain;
29+
1. The inputs are decorrelated with their sample mean and covariance (and projected to low dimensional subspace if necessary) i.e PCA
30+
2. The scaled inputs are then subject to a "Robust" univariate scaling, mapping 1st-3rd quartiles to [0,1]
31+
3. The outputs are decorrelated using an "output structure matrix" (provided to the emulator `output_structure_matrix=`). Furthermore, apply a dimension-reduction to a space that retains 95% of the total variance.
32+
4. In the reduced input-output space, a canonical correlation analysis is performed. Data is oriented and reduced (if necessary) maximize the joint correlation between inputs and outputs.
33+
34+
!!! note "Default Encoder schedule"
35+
The current default encoder schedule applies `decorrelate_structure_mat()` if a structure matrix (input or output) is provided, else it applies `decorrelate_sample_cov()`.
36+
37+
!!! note "Switch encoding off"
38+
To ensure that no encoding is happening, the user must pass in an empty schedule `encoder_schedule = []`
39+
40+
41+
## Creating an emulator with a schedule
42+
43+
The schedule is then passed into the Emulator, along with the data and desired structure matrices
44+
```julia
45+
emulator = Emulator(
46+
machine_learning_tool,
47+
input_output_pairs;
48+
output_structure_matrix = obs_noise_cov,
49+
encoder_schedule = complex_schedule,
50+
)
51+
```
52+
Note that due to the item `(decorrelate_structure_mat(retain_var=0.95), "out")` in the schedule, we must provide the `output_structure_matrix`.
53+
54+
# Types of data processors
55+
56+
We currently provide two main types of data processing: the `DataContainerProcessor` and `PairedDataContainerProcessor`.
57+
58+
The `DataContainerProcessor` encodes "input" data agnostic of the "output" data, and vice versa, examples of current implementations are:
59+
- `UnivariateAffineScaling`: such as `quartile_scale()`, `minmax_scale()`, and `zscore_scale()`, which apply some basic univariate scaling to the data in each dimension
60+
- `Decorrelator`: such as `decorrelate_structure_mat()` and `decorrelate_sample_cov()`, or `decorrelate()` which perform [(truncated-)PCA](https://en.wikipedia.org/wiki/Singular_value_decomposition) using either the sample-estimated or user-provided covariance matrices (or their sum).
61+
62+
The `PairedDataContainerProcessor` encodes inputs (or outputs) using information of the both inputs and outputs in pairs
63+
- `CanonicalCorrelation` - constructed with `canonical_correlation()`, which performs [canonical correlation analysis](https://en.wikipedia.org/wiki/Canonical_correlation) to process the pairs. In effect this performs PCA on the cross-correlation from input and output samples.
64+
- [Coming soon] `LikelihoodInformed` - this will use the data or likelihood from the inverse problem at hand to build diagnostic matrices that are used to find informative directions for dimension reduction (e.g., [Cui, Zahm 2021](http://doi.org/10.1088/1361-6420/abeafb), [Baptista, Marzouk, Zahm 2022](https://arxiv.org/abs/2207.08670))
65+
66+
This is an extensible framework, and so new data processors can be added to this library.
67+
68+
!!! note "Some experiments"
69+
Some effects of the data processing (with older API) are outlined in a practical setting in the results and appendices of [Howland, Dunbar, Schneider, (2022)](https://doi.org/10.1029/2021MS002735).

docs/src/emulate.md

Lines changed: 3 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -17,14 +17,11 @@ Wrapping a predefined machine learning tool, e.g. a Gaussian process `gauss_proc
1717
emulator = Emulator(
1818
gauss_proc,
1919
input_output_pairs; # optional arguments after this
20-
obs_noise_cov = Γy,
21-
normalize_inputs = true,
22-
standardize_outputs = true,
23-
standardize_outputs_factors = factor_vector,
24-
retained_svd_frac = 0.95,
20+
output_structure_matrix = Γy,
21+
encoder_schedule = encoder_schedule,
2522
)
2623
```
27-
The optional arguments above relate to the data processing.
24+
The optional arguments above relate to the data processing, which is described [here](@ref data-proc)
2825

2926
### Emulator Training
3027

@@ -40,34 +37,6 @@ y, cov = Emulator.predict(emulator, new_inputs)
4037
```
4138
This returns both a mean value and a covariance.
4239

43-
44-
## Data processing
45-
46-
Some effects of the following are outlined in a practical setting in the results and appendices of [Howland, Dunbar, Schneider, (2022)](https://doi.org/10.1029/2021MS002735).
47-
48-
### Diagonalization and output dimension reduction
49-
50-
This arises from the optional arguments
51-
- `obs_noise_cov = Γy` (default: `nothing`)
52-
We always use singular value decomposition to diagonalize the output space, requiring output covariance `Γy`. *Why?* If we need to train a $$\mathbb{R}^{10} \to \mathbb{R}^{100}$$ emulator, diagonalization allows us to instead train 100 $$\mathbb{R}^{10} \to \mathbb{R}^{1}$$ emulators (far cheaper).
53-
- `retained_svd_frac = 0.95` (default `1.0`)
54-
Performance is increased further by throwing away less informative output dimensions, if 95% of the information (i.e., variance) is in the first 40 diagonalized output dimensions then setting `retained_svd_frac=0.95` will train only 40 emulators.
55-
56-
!!! note
57-
Diagonalization is an approximation. It is however a good approximation when the observational covariance varies slowly in the parameter space.
58-
!!! warn
59-
Severe approximation errors can occur if `obs_noise_cov` is not provided.
60-
61-
62-
### Normalization and standardization
63-
64-
This arises from the optional arguments
65-
- `normalize_inputs = true` (default: `true`)
66-
We normalize the input data in a standard way by centering, and scaling with the empirical covariance
67-
- `standardize_outputs = true` (default: `false`)
68-
- `standardize_outputs_factors = factor_vector` (default: `nothing`)
69-
To help with poor conditioning of the covariance matrix, users can also standardize each output dimension with by a multiplicative factor given by the elements of `factor_vector`.
70-
7140
## [Modular interface](@id modular-interface)
7241

7342
Developers may contribute new tools by performing the following

examples/Cloudy/Cloudy_calibrate.jl

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -160,13 +160,19 @@ dummy = ones(n_params)
160160
dist_type = ParticleDistributions.GammaPrimitiveParticleDistribution(dummy...)
161161
model_settings = DynamicalModel.ModelSettings(kernel, dist_type, moments, tspan)
162162
# EKI iterations
163+
n_iter = N_iter
163164
for n in 1:N_iter
164165
# Return transformed parameters in physical/constrained space
165166
ϕ_n = get_ϕ_final(priors, ekiobj)
166167
# Evaluate forward map
167168
G_n = [DynamicalModel.run_dyn_model(ϕ_n[:, i], model_settings) for i in 1:N_ens]
168169
G_ens = hcat(G_n...) # reformat
169-
EnsembleKalmanProcesses.update_ensemble!(ekiobj, G_ens)
170+
terminate = EnsembleKalmanProcesses.update_ensemble!(ekiobj, G_ens)
171+
if !isnothing(terminate)
172+
n_iter = n - 1
173+
break
174+
end
175+
170176
end
171177

172178

@@ -202,7 +208,7 @@ save(
202208
gr(size = (1200, 400))
203209

204210
u_init = get_u_prior(ekiobj)
205-
anim_eki_unconst_cloudy = @animate for i in 1:(N_iter - 1)
211+
anim_eki_unconst_cloudy = @animate for i in 1:(n_iter - 1)
206212
u_i = get_u(ekiobj, i)
207213

208214
p1 = plot(u_i[1, :], u_i[2, :], seriestype = :scatter, xlims = extrema(u_init[1, :]), ylims = extrema(u_init[2, :]))
@@ -259,7 +265,7 @@ gif(anim_eki_unconst_cloudy, joinpath(output_directory, "cloudy_eki_unconstr.gif
259265

260266
# Plots in the constrained space
261267
ϕ_init = transform_unconstrained_to_constrained(priors, u_init)
262-
anim_eki_cloudy = @animate for i in 1:(N_iter - 1)
268+
anim_eki_cloudy = @animate for i in 1:(n_iter - 1)
263269
ϕ_i = get_ϕ(priors, ekiobj, i)
264270

265271
p1 = plot(ϕ_i[1, :], ϕ_i[2, :], seriestype = :scatter, xlims = extrema(ϕ_init[1, :]), ylims = extrema(ϕ_init[2, :]))

examples/Cloudy/Cloudy_emulate_sample.jl

Lines changed: 21 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -16,16 +16,7 @@ using CairoMakie, PairPlots
1616
using CalibrateEmulateSample.Emulators
1717
using CalibrateEmulateSample.MarkovChainMonteCarlo
1818
using CalibrateEmulateSample.Utilities
19-
using EnsembleKalmanProcesses
20-
using EnsembleKalmanProcesses.ParameterDistributions
21-
using EnsembleKalmanProcesses.DataContainers
22-
23-
function get_standardizing_factors(data::Array{FT, 2}) where {FT}
24-
# Input: data size: N_data x N_ensembles
25-
# Ensemble median of the data
26-
norm_factor = median(data, dims = 2) # N_data x 1 array
27-
return norm_factor
28-
end
19+
using CalibrateEmulateSample.EnsembleKalmanProcesses
2920

3021
################################################################################
3122
# #
@@ -100,7 +91,7 @@ function main()
10091
cases = [
10192
"rf-scalar",
10293
"gp-gpjl", # Veeeery slow predictions
103-
"rf-nosvd-nonsep",
94+
"rf-nonsep",
10495
]
10596

10697
# Specify cases to run (e.g., case_mask = [2] only runs the second case)
@@ -116,13 +107,13 @@ function main()
116107
"verbose" => true,
117108
"scheduler" => DataMisfitController(terminate_at = 100.0),
118109
"cov_sample_multiplier" => 1.0,
119-
"n_iteration" => 20,
110+
"n_iteration" => 10,
120111
)
121112

122113
# We use the same input-output-pairs and normalization factors for
123114
# Gaussian Process and Random Feature cases
124115
input_output_pairs = get_training_points(ekiobj, length(get_u(ekiobj)) - 2)
125-
norm_factors = get_standardizing_factors(get_outputs(input_output_pairs))
116+
test_pairs = get_training_points(ekiobj, (length(get_u(ekiobj)) - 1):(length(get_u(ekiobj)) - 1))
126117
for case in cases[case_mask]
127118

128119
println(" ")
@@ -140,12 +131,9 @@ function main()
140131
# Define machine learning tool
141132
mlt = GaussianProcess(gppackage; kernel = gp_kernel, prediction_type = pred_type, noise_learn = false)
142133

143-
decorrelate = true
144-
standardize_outputs = true
145-
146134
elseif case == "rf-scalar"
147135

148-
kernel_rank = 3
136+
kernel_rank = 1
149137
kernel_structure = SeparableKernel(LowRankFactor(kernel_rank, nugget), OneDimFactor())
150138

151139
# Define machine learning tool
@@ -156,10 +144,7 @@ function main()
156144
optimizer_options = optimizer_options,
157145
)
158146

159-
decorrelate = true
160-
standardize_outputs = true
161-
162-
elseif case == "rf-nosvd-nonsep"
147+
elseif case == "rf-nonsep"
163148

164149
# Define machine learning tool
165150
kernel_rank = 4
@@ -171,34 +156,32 @@ function main()
171156
optimizer_options = optimizer_options,
172157
)
173158

174-
# Vector RF does not require decorrelation of outputs
175-
decorrelate = false
176-
standardize_outputs = false
177-
178-
179159
else
180160
error("Case $case is not implemented yet.")
181161

182162
end
183163

184-
# The data processing normalizes input data, and decorrelates
185-
# output data with information from Γy, if required
186-
# Note: The `standardize_outputs_factors` are only used under the
187-
# condition that `standardize_outputs` is true.
164+
# Data processing
165+
if case == "rf-nonsep"
166+
encoder_schedule = []
167+
else
168+
encoder_schedule = (decorrelate_structure_mat(), "in_and_out")
169+
end
170+
171+
# build emulator
188172
emulator = Emulator(
189173
mlt,
190-
input_output_pairs,
191-
obs_noise_cov = Γy,
192-
decorrelate = decorrelate,
193-
standardize_outputs = standardize_outputs,
194-
standardize_outputs_factors = vcat(norm_factors...),
174+
input_output_pairs;
175+
input_structure_matrix = cov(priors),
176+
output_structure_matrix = Γy,
177+
encoder_schedule = encoder_schedule,
195178
)
196179

197180
optimize_hyperparameters!(emulator)
198181

199182
# Check how well the emulator predicts on the true parameters
200183
y_mean, y_var = Emulators.predict(emulator, reshape(θ_true, :, 1); transform_to_real = true)
201-
184+
y_mean_test, y_var_test = Emulators.predict(emulator, get_inputs(test_pairs); transform_to_real = true)
202185
println("Emulator ($(case)) prediction on true parameters: ")
203186
println(vec(y_mean))
204187
println("true data: ")
@@ -207,6 +190,8 @@ function main()
207190
println(sqrt.(diag(y_var[1], 0)))
208191
println("Emulator ($(case)) MSE (truth): ")
209192
println(mean((truth_sample - vec(y_mean)) .^ 2))
193+
println("Emulator ($(case)) MSE (next ensemble): ")
194+
println(mean((get_outputs(test_pairs) - y_mean_test) .^ 2))
210195

211196

212197
###

examples/Cloudy/Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,4 +17,4 @@ StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
1717

1818
[compat]
1919
Cloudy = "0.2"
20-
julia = "~1.6"
20+
julia = "1.6"

0 commit comments

Comments
 (0)