diff --git a/README.md b/README.md index 7796296e1..d8ecee242 100644 --- a/README.md +++ b/README.md @@ -36,3 +36,12 @@ Implements a derivative-free machine-learning-accelerated pipeline for uncertain [dlt-url]: https://juliapkgstats.com/pkg/CalibrateEmulateSample.c +### Quick links + +- [How do I build prior distributions?](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/parameter_distributions/) +- [How do I build good observational noise covariances](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/observations/) +- [What ensemble size should I take? Which process should I use? What is the recommended configuration?](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/defaults/) +- [Where can I walk through the simple example?](https://clima.github.io/CalibrateEmulateSample.jl/dev/examples/sinusoid_example/) +- [What is the `EnsembleKalmanProcesses.jl` package?](https://clima.github.io/CalibrateEmulateSample.jl/dev/calibrate/) +- [What are the recommendations/defaults for dimension reduction or data processing?](https://clima.github.io/CalibrateEmulateSample.jl/dev/data_processing/) +- [How to I plot or interpret the posterior distribution?](https://clima.github.io/CalibrateEmulateSample.jl/dev/sample/) diff --git a/docs/src/GaussianProcessEmulator.md b/docs/src/GaussianProcessEmulator.md index cbb5e8b22..ead5a9642 100644 --- a/docs/src/GaussianProcessEmulator.md +++ b/docs/src/GaussianProcessEmulator.md @@ -103,6 +103,45 @@ gauss_proc = GaussianProcess( ``` You can also combine multiple ScikitLearn kernels via linear operations in the same way as above. +## [AGPJL](@id agpjl) + +Autodifferentiable emulators are used by our differentiable samplers. Currently the only support for autodifferentiable Gaussian process emulators in Julia (within the `predict()` method, not hyperparameter optimization) is to use `AbstractGPs.jl`. As `AbstractGPs.jl` has no optimization routines for kernels, we instead apply the following (temporary) recipe: +- Create and optimize a `GPJL` emulator and default kernel. +```julia +gp_jl = GaussianProcess(GPJL(); noise_learn=true, gpjl_kwargs...) # must! use default kernel +em = Emulator(gp_jl, iopairs) +optimize_hyperparameters!(em) # updates gp_jl +``` +- Create the Kernel parameters as a vect-of-dict with +```julia +kernel_params = [ + Dict( + "log_rbf_len" => model_params[1:end-2] # input-dim Vector, + "log_std_sqexp" => model_params[end-1] # Float, + "log_std_noise" => model_params[end],# Float, + ) +for model_params in get_params(gp_jl)] +``` +- Note: `get_params(gp_jl)` returns `output_dim`-vector where each entry is [a, b, c] with: + - a is the `rbf_len`: lengthscale parameters for SEArd kernel [input_dim]- length Vector + - b is the `log_std_sqexp` of the SQexp kernel Float + - c is the `log_std_noise` of the noise kernel Float +- Build a new gaussian process with `AGPJL`, use same keyword arguments as with `GPJL()` and add the kernel_params kwarg into the emulator +```julia +agp_jl = GaussianProcess(AGPJL(); noise_learn=true, gpjl_kwargs...) +em = Emulator(agp_jl, iopairs; kernel_params=kernel_params) +# no call to optimize_hyperparameters +``` + +We would be keen to see contributions to our codebase to improve this interface (or to perform hyperparameter optimization with AGP directly). + +!!! note "Other notes about the recipe" + 1. As stated, this only works for the default kernel option. + 2. If `noise_learn=false` then `log_std_noise` does not appear, the `model_params` indexing to the other hyperparameters should be adjusted + 2. For one-dimensional output, use `opt_params = Emulators.get_params(gauss_proc)[1]` to get the parameters + + + # Learning additional white noise Often it is useful to learn the discrepancy between the Gaussian process prediction and the data, by learning additional white noise. Though one often knows, and provides, the discrepancy between the true model and data with an observational noise covariance; the additional white kernel can help account for approximation error from the selected Gaussian process kernel and the true model. This is added with the Boolean keyword `noise_learn` when initializing the Gaussian process. The default is true. diff --git a/docs/src/assets/GFunction_sens_RF-scalar_10.png b/docs/src/assets/GFunction_sens_RF-scalar_10.png index 2a7adb40f..5205fe7ff 100644 Binary files a/docs/src/assets/GFunction_sens_RF-scalar_10.png and b/docs/src/assets/GFunction_sens_RF-scalar_10.png differ diff --git a/docs/src/assets/GFunction_sens_RF-scalar_3.png b/docs/src/assets/GFunction_sens_RF-scalar_3.png index 754280d1c..8b9948037 100644 Binary files a/docs/src/assets/GFunction_sens_RF-scalar_3.png and b/docs/src/assets/GFunction_sens_RF-scalar_3.png differ diff --git a/docs/src/assets/GFunction_slices_RF-scalar_10.png b/docs/src/assets/GFunction_slices_RF-scalar_10.png index 08d0bd81e..803a117ad 100644 Binary files a/docs/src/assets/GFunction_slices_RF-scalar_10.png and b/docs/src/assets/GFunction_slices_RF-scalar_10.png differ diff --git a/docs/src/assets/GFunction_slices_RF-scalar_3.png b/docs/src/assets/GFunction_slices_RF-scalar_3.png index 727f2a594..dbdf7be02 100644 Binary files a/docs/src/assets/GFunction_slices_RF-scalar_3.png and b/docs/src/assets/GFunction_slices_RF-scalar_3.png differ diff --git a/docs/src/assets/GP_l63_attr.png b/docs/src/assets/GP_l63_attr.png index eea033d22..afb762672 100644 Binary files a/docs/src/assets/GP_l63_attr.png and b/docs/src/assets/GP_l63_attr.png differ diff --git a/docs/src/assets/GP_l63_cdfs.png b/docs/src/assets/GP_l63_cdfs.png index b5f98a1ed..60ed4773a 100644 Binary files a/docs/src/assets/GP_l63_cdfs.png and b/docs/src/assets/GP_l63_cdfs.png differ diff --git a/docs/src/assets/GP_l63_pdf.png b/docs/src/assets/GP_l63_pdf.png index 65427873d..5efb50079 100644 Binary files a/docs/src/assets/GP_l63_pdf.png and b/docs/src/assets/GP_l63_pdf.png differ diff --git a/docs/src/assets/GP_l63_test.png b/docs/src/assets/GP_l63_test.png index 254d559f8..7693721d3 100644 Binary files a/docs/src/assets/GP_l63_test.png and b/docs/src/assets/GP_l63_test.png differ diff --git a/docs/src/assets/RF-nonsep_l63_attr.png b/docs/src/assets/RF-nonsep_l63_attr.png new file mode 100644 index 000000000..b319a1e36 Binary files /dev/null and b/docs/src/assets/RF-nonsep_l63_attr.png differ diff --git a/docs/src/assets/RF-nonsep_l63_cdfs.png b/docs/src/assets/RF-nonsep_l63_cdfs.png new file mode 100644 index 000000000..2921cb133 Binary files /dev/null and b/docs/src/assets/RF-nonsep_l63_cdfs.png differ diff --git a/docs/src/assets/RF-nonsep_l63_pdf.png b/docs/src/assets/RF-nonsep_l63_pdf.png new file mode 100644 index 000000000..1a3b8e66a Binary files /dev/null and b/docs/src/assets/RF-nonsep_l63_pdf.png differ diff --git a/docs/src/assets/RF-nonsep_l63_test.png b/docs/src/assets/RF-nonsep_l63_test.png new file mode 100644 index 000000000..43ae307ca Binary files /dev/null and b/docs/src/assets/RF-nonsep_l63_test.png differ diff --git a/docs/src/assets/RF-nosvd-nonsep_l63_attr.png b/docs/src/assets/RF-nosvd-nonsep_l63_attr.png deleted file mode 100644 index 628489cfe..000000000 Binary files a/docs/src/assets/RF-nosvd-nonsep_l63_attr.png and /dev/null differ diff --git a/docs/src/assets/RF-nosvd-nonsep_l63_cdfs.png b/docs/src/assets/RF-nosvd-nonsep_l63_cdfs.png deleted file mode 100644 index e890d2fc7..000000000 Binary files a/docs/src/assets/RF-nosvd-nonsep_l63_cdfs.png and /dev/null differ diff --git a/docs/src/assets/RF-nosvd-nonsep_l63_pdf.png b/docs/src/assets/RF-nosvd-nonsep_l63_pdf.png deleted file mode 100644 index d42dc47e1..000000000 Binary files a/docs/src/assets/RF-nosvd-nonsep_l63_pdf.png and /dev/null differ diff --git a/docs/src/assets/RF-nosvd-nonsep_l63_test.png b/docs/src/assets/RF-nosvd-nonsep_l63_test.png deleted file mode 100644 index 7a4ffc0fa..000000000 Binary files a/docs/src/assets/RF-nosvd-nonsep_l63_test.png and /dev/null differ diff --git a/docs/src/assets/sinusoid_GP_emulator_contours.png b/docs/src/assets/sinusoid_GP_emulator_contours.png index f16d4e083..057885b84 100644 Binary files a/docs/src/assets/sinusoid_GP_emulator_contours.png and b/docs/src/assets/sinusoid_GP_emulator_contours.png differ diff --git a/docs/src/assets/sinusoid_GP_errors_contours.png b/docs/src/assets/sinusoid_GP_errors_contours.png index 8cf5a8211..7a273052e 100644 Binary files a/docs/src/assets/sinusoid_GP_errors_contours.png and b/docs/src/assets/sinusoid_GP_errors_contours.png differ diff --git a/docs/src/assets/sinusoid_MCMC_hist_GP.png b/docs/src/assets/sinusoid_MCMC_hist_GP.png index 3f87b4827..540c710a7 100644 Binary files a/docs/src/assets/sinusoid_MCMC_hist_GP.png and b/docs/src/assets/sinusoid_MCMC_hist_GP.png differ diff --git a/docs/src/assets/sinusoid_MCMC_hist_RF.png b/docs/src/assets/sinusoid_MCMC_hist_RF.png index 0aa9db9fd..555d4fdc0 100644 Binary files a/docs/src/assets/sinusoid_MCMC_hist_RF.png and b/docs/src/assets/sinusoid_MCMC_hist_RF.png differ diff --git a/docs/src/assets/sinusoid_RF_emulator_contours.png b/docs/src/assets/sinusoid_RF_emulator_contours.png index 04b484021..4f66f0fc1 100644 Binary files a/docs/src/assets/sinusoid_RF_emulator_contours.png and b/docs/src/assets/sinusoid_RF_emulator_contours.png differ diff --git a/docs/src/assets/sinusoid_RF_errors_contours.png b/docs/src/assets/sinusoid_RF_errors_contours.png index b864e1d78..a0730f508 100644 Binary files a/docs/src/assets/sinusoid_RF_errors_contours.png and b/docs/src/assets/sinusoid_RF_errors_contours.png differ diff --git a/docs/src/assets/sinusoid_eki_pairs.png b/docs/src/assets/sinusoid_eki_pairs.png index 451d31996..8ffdfd73f 100644 Binary files a/docs/src/assets/sinusoid_eki_pairs.png and b/docs/src/assets/sinusoid_eki_pairs.png differ diff --git a/docs/src/assets/sinusoid_groundtruth_contours.png b/docs/src/assets/sinusoid_groundtruth_contours.png index 289044581..e30c34dde 100644 Binary files a/docs/src/assets/sinusoid_groundtruth_contours.png and b/docs/src/assets/sinusoid_groundtruth_contours.png differ diff --git a/docs/src/assets/sinusoid_posterior_GP.png b/docs/src/assets/sinusoid_posterior_GP.png index 046765e89..4ae138967 100644 Binary files a/docs/src/assets/sinusoid_posterior_GP.png and b/docs/src/assets/sinusoid_posterior_GP.png differ diff --git a/docs/src/assets/sinusoid_posterior_RF.png b/docs/src/assets/sinusoid_posterior_RF.png index 2f022118e..f07f9bae8 100644 Binary files a/docs/src/assets/sinusoid_posterior_RF.png and b/docs/src/assets/sinusoid_posterior_RF.png differ diff --git a/docs/src/assets/sinusoid_prior.png b/docs/src/assets/sinusoid_prior.png index ae7e41d1f..42c1b790c 100644 Binary files a/docs/src/assets/sinusoid_prior.png and b/docs/src/assets/sinusoid_prior.png differ diff --git a/docs/src/assets/sinusoid_true_vs_observed_signal.png b/docs/src/assets/sinusoid_true_vs_observed_signal.png index c01062fda..f1ff25658 100644 Binary files a/docs/src/assets/sinusoid_true_vs_observed_signal.png and b/docs/src/assets/sinusoid_true_vs_observed_signal.png differ diff --git a/docs/src/calibrate.md b/docs/src/calibrate.md index 7874430cc..44fad7bf1 100644 --- a/docs/src/calibrate.md +++ b/docs/src/calibrate.md @@ -1,4 +1,4 @@ -# The Calibrate stage +# [The Calibrate stage](@id calibrate) Calibration of the computer model entails finding an optimal parameter ``\theta^*`` that maximizes the posterior probability diff --git a/docs/src/data_processing.md b/docs/src/data_processing.md index 1e17464dc..e1db4e328 100644 --- a/docs/src/data_processing.md +++ b/docs/src/data_processing.md @@ -1,14 +1,81 @@ # [Data Processing and Dimension Reduction](@id data-proc) ## Overview -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 +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. To make results as transparent and reproducible as possible we provide a flexible (and extensible!) framework to create encoders and decoders for this purpose. + +The framework works as follows - An `encoder_schedule` defines the type of processing to be applied to input and/or output spaces -- The `encoder_schedule` is passed into the `Emulator` where it is initialized and stored. -- The encoder will be used automatically to encode training data and predictions within the Emulate and Sample routines. +- The `encoder_schedule` is passed into the `Emulator` where it is initialized and stored. Sometimes it will require additional information, passed in as `encoder_kwargs`. +- The encoder will be used automatically to encode training data, covariance matrices, and predictions within the Emulate and Sample routines. An external API is also available using the `encode_data`, and `encode_structure_matrix` methods if needed. -## Define an encoder schedule +## Defaults and recommendations +### Default schedule (no explicit dimension reduction) +When no schedule is provided, i.e. +```julia +emulator = Emulator(machine_learning_tool, input_output_pairs) +``` +The default schedule under-the-hood, is given by +```julia +schedule = (decorrelate_sample_cov(), "in_and_out") +``` +If the user only provides the observational noise covariance, +```julia +emulator = Emulator( + machine_learning_tool, + input_output_pairs; + encoder_kwargs = (; obs_noise_cov = obs_noise_cov), +) +``` +The default schedule under-the-hood, is given by +```julia +schedule = [ + (decorrelate_sample_cov(), "in"), + (decorrelate_structure_mat(), "out"), # uses obs_noise_cov +] +``` + +!!! note "Switch encoding off" + To ensure that no encoding is happening, the user must pass in an empty schedule `encoder_schedule = []` + +### Recommended: schedule for PCA dimension reduction + +Our recommendeded family to balance efficiency and reduce dimension is to use the `retain_var` kwargs. +```julia +retain_var_in = 0.99 # reduce dimension retaining 99% of input variance +retain_var_out = 0.95 # reduce dimension retaining 95% of output variance +encoder_schedule = [ + (decorrelate_sample_cov(retain_var = retain_var_in), "in"), + (decorrelate_structure_mat(retain_var = retain_var_out), "out"), +] +encoder_kwargs = (; obs_noise_cov = obs_noise_cov) + +emulator = Emulator( + machine_learning_tool, + input_output_pairs; + encoder_schedule = encoder_schedule, + encoder_kwargs = encoder_kwargs, +) +``` +## Recommended: Get `encoder_kwargs` from `ObservationSeries` etc. + +To transition more smoothly from the `EnsembleKalmanProcesses` infrastructure, one can also get the kwargs in the desired format from any [`Observation` and `ObservationSeries`](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/observations/) objects, as well as [`ParameterDistribution`](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/parameter_distributions/) objects used in the `EnsembleKalmanProcesses.jl` package. +```julia +# `prior::ParameterDistribution` +input_kwargs = get_kwargs_from(prior) + +# `observation_series::ObservationSeries` +output_kwargs = get_kwargs_from(observation_series) + +encoder_kwargs = merge(input_kwargs, output_kwargs) +``` +!!! note "Observational noise notes" + 1. The `obs_noise_cov` object does not need to be a constructed `AbstractMatrix`, rather it can be any type of compactly stored matrix compatible with `EnsembleKalmanProcesses` (e.g., [here](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/observations/)). + 2. When `ObservationSeries` contains multiple observations, it is assumed that the covariance of the noise of all observations are the same. (If they are not, only the first will be taken by default in `get_kwargs_from(...)` + + +## Building and interpreting encoder schedules 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] ```julia @@ -34,10 +101,6 @@ In this (rather unrealistic) chain; !!! note "Default Encoder schedule" The current default encoder schedule applies `decorrelate_structure_mat()` if a structure matrix (input or output) is provided, else it applies `decorrelate_sample_cov()`. -!!! note "Switch encoding off" - To ensure that no encoding is happening, the user must pass in an empty schedule `encoder_schedule = []` - - ## Creating an emulator with a schedule The schedule is then passed into the Emulator, along with the data and desired structure matrices @@ -51,6 +114,26 @@ emulator = Emulator( ``` Note that due to the item `(decorrelate_structure_mat(retain_var=0.95), "out")` in the schedule, we must provide an output structure matrix. In this case, we provide `obs_noise_cov`. +### Additional API +Once created, the user does not need to worry about encoding/decoding. However, we do provide an external API that can be used to encode both new data vectors (or matrix-of-columns), and structure matrices (i.e. covariances) for validation or error diagnosis. +```julia +# encoding new data: Vector, `DataContainer`, or Matrix viewed as columns +enc_in_data = encode_data(emulator, new_input_data, "in") +enc_out_data = encode_data(emulator, new_output_data, "out") + +# decoding it +dec_in_data = decode_data(emulator, enc_in_data, "in") +dec_out_data = decode_data(emulator, enc_out_data, "out") + +# new structure matrix: Matrix (typically acting as a covariance) +enc_in_cov = encode_structure_matrix(emulator, new_input_covariance, "in") +enc_out_cov = encode_structure_matrix(emulator, new_output_covariance, "out") + +dec_in_cov = decode_structure_matrix(emulator, enc_in_cov, "in") +dec_out_cov = decode_structure_matrix(emulator, enc_out_cov, "out") +``` + + # Types of data processors We currently provide two main types of data processing: the `DataContainerProcessor` and `PairedDataContainerProcessor`. @@ -61,9 +144,17 @@ The `DataContainerProcessor` encodes "input" data agnostic of the "output" data, The `PairedDataContainerProcessor` encodes inputs (or outputs) using information of the both inputs and outputs in pairs - `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. -- [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)) +- [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. In particular we build generalizations of current frameworks (e.g., [Cui, Zahm 2021](http://doi.org/10.1088/1361-6420/abeafb), [Baptista, Marzouk, Zahm 2022](https://arxiv.org/abs/2207.08670)) to use the latest EKP iterations. This is an extensible framework, and so new data processors can be added to this library. !!! note "Some experiments" 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). + +## A note on `LinearMaps` and the internals for encoding large observations + +To handle the heterogeneity and possible massive size of matrices involved with encoding observations, all `obs_noise_cov` and encoder-decoder objects are converted into linear maps, using [`LinearMaps.jl`](https://julialinearalgebra.github.io/LinearMaps.jl/stable/). For high-dimensional inputs or outputs (>3000), these are used in within iterative or (stochastically) approximate matrix-free methods, thus giving some possible additional error to gain better scaling in dimension. Such balances are controllable by the user with keywords that can be passed to certain encoders in the schedule. + + + + diff --git a/docs/src/examples/Cloudy_example.md b/docs/src/examples/Cloudy_example.md index f34f2e51b..0c227f57d 100755 --- a/docs/src/examples/Cloudy_example.md +++ b/docs/src/examples/Cloudy_example.md @@ -35,7 +35,7 @@ The input to the CES algorithm consists of data ``y``, the observational covaria CES to learn model parameters. -### Running the example +# Running the example `Cloudy_calibrate.jl` performs the calibration using ensemble Kalman inversion; `Cloudy_emulate_sample.jl` fits an emulator and uses it to sample the @@ -46,11 +46,11 @@ include("Cloudy_calibrate.jl") include("Cloudy_emulate_sample.jl") ``` -### Walkthrough of the code: `Cloudy_calibrate.jl` +# Walkthrough of the code: `Cloudy_calibrate.jl` This file performs the calibration stage of CES. -#### Import packagages and modules +## Import packagages and modules First we load standard packages, @@ -92,7 +92,7 @@ that runs Cloudy with a given instance of the parameter vector we want to learn. include("DynamicalModel.jl") ``` -#### Define the true parameters +### Define the true parameters We define the true parameters---they are known here because this is a known truth example. Knowing the true parameters will allow us to assess how well @@ -113,7 +113,7 @@ k_true = 0.0817 # shape parameter of Gamma distribution dist_true = ParticleDistributions.GammaPrimitiveParticleDistribution(ϕ_true...) ``` -#### Define priors for the parameters +### Define priors for the parameters As we are working with Bayesian methods, we treat the parameters we want to learn as random variables whose prior distributions we specify here. The prior @@ -141,7 +141,7 @@ p = plot(priors, constrained=false) ![priors](../assets/cloudy_priors.png) -#### Generate (synthetic) observations +### Generate (synthetic) observations We generate synthetic observations by running Cloudy 100 times with the true parameters (i.e., with the true initial Gamma distribution of droplet masses) and then adding noise to simulate measurement error. @@ -175,7 +175,7 @@ truth = Observation( ) ``` -#### Perform ensemble Kalman inversion +### Perform ensemble Kalman inversion We sample the initial ensemble from the prior and create the `EnsembleKalmanProcess` object as an ensemble Kalman inversion (EKI) algorithm @@ -216,7 +216,7 @@ for n in 1:N_iter end ``` -#### Visualize and store the results of the calibration +### Visualize and store the results of the calibration The optimal parameter vector determined by the ensemble Kalman inversion is the ensemble mean of the particles after the last iteration, which is printed to @@ -233,12 +233,12 @@ computational (unconstrained) and physical (constrained) spaces. ![eki_iterations_animation](../assets/cloudy_eki.gif) -### Walkthrough of the code: `Cloudy_emulate_sample.jl` +## Walkthrough of the code: `Cloudy_emulate_sample.jl` This file performs the emulation and sampling stages of the CES algorithm. -#### Import packages and modules +### Import packages and modules First, we import some standard packages ```julia @@ -264,7 +264,7 @@ using EnsembleKalmanProcesses.ParameterDistributions using EnsembleKalmanProcesses.DataContainers ``` -#### Load the calibration results +### Load the calibration results We will train the emulator on the input-output pairs we obtained during the calibration. They are stored within the `EnsembleKalmanProcess` object @@ -381,34 +381,25 @@ standardize_outputs = false We construct the emulator using the input-output pairs obtained in the calibration stage (note that we're not using all available input-output pairs---using all of them may not give the best results, especially if the EKI -parameter converges rapidly and then "stays in the same place" during the remaining iterations). For the `gp-gpjl` and `rf-scalar` cases, we want the output -data to be decorrelated with information from Γy, but for the vector RF case -decorrelation is not required. +parameter converges rapidly and then "stays in the same place" during the remaining iterations). For the `gp-gpjl` and `rf-scalar` cases, create the data processor and build the emulator as follows ``` input_output_pairs = get_training_points(ekiobj, length(get_u(ekiobj))-2) -# Use the medians of the outputs as standardizing factors -norm_factors = get_standardizing_factors( - get_outputs(input_output_pairs) -) +# The data processing decorrelates inputs and outputs using the prior and noise covariances respectively +encoder_schedule = (decorrelate_structure_mat(), "in_and_out") +encoder_kwargs = (; prior_cov = cov(priors), obs_noise_cov = Γy) -# The data processing normalizes input data, and decorrelates -# output data with information from Γy, if required -# Note: The `standardize_outputs_factors` are only used under the -# condition that `standardize_outputs` is true. emulator = Emulator( mlt, - input_output_pairs, - decorrelate = decorrelate, - obs_noise_cov = Γy, - standardize_outputs = true, - standardize_outputs_factors = vcat(norm_factors...), + input_output_pairs; + encoder_schedule = encoder_schedule, + encoder_kwargs = encoder_kwargs, ) ``` -#### Train the emulator +## Train the emulator The emulator is trained when we combine the machine learning tool and the data into the `Emulator` above. We must also optimize the hyperparameters: @@ -438,7 +429,7 @@ println(truth_sample) # what was used as truth The emulator predicts both a mean value and a covariance. -### Sample the posterior distributions of the parameters +## Sample the posterior distributions of the parameters The last step is to plug the emulator into an MCMC algorithm, which is then used to produce samples from the posterior distribution of the parameters. Essentially, the emulator acts as a stand-in for the original forward model (which in most cases of interest is computationally expensive to run) during the MCMC sampling process. @@ -495,7 +486,7 @@ the peak of the posterior distribution should be located near the true values, indicating a high-quality estimation. Additionally, visualizing the prior distributions alongside the posteriors shows the distributional change effected by the Bayesian learning process. -### Results +## Results We first produce pair plots (also known as corner plots or scatter plot matrices) to visualize the posterior parameter distributions as a grid of histograms. Recall that the task was to solve the inverse problem of finding the parameters ``N_{0, 0}``, ``k_0``, and ``\theta_0``, which define a gamma distribution of droplet masses in Cloudy at time ``t = 0``. diff --git a/docs/src/examples/emulators/global_sens_analysis.md b/docs/src/examples/emulators/global_sens_analysis.md index 79f22a77e..3bb794c16 100644 --- a/docs/src/examples/emulators/global_sens_analysis.md +++ b/docs/src/examples/emulators/global_sens_analysis.md @@ -117,9 +117,11 @@ mlt = GaussianProcess( noise_learn = false, ) ``` -We finish by building the emulator object +We finish by defining the data processing, and building the emulator object ```julia -emulator = Emulator(mlt, iopairs; obs_noise_cov = Γ * I, decorrelate = decorrelate) +encoder_schedule = (decorrelate_structure_mat(), "out") +encoder_kwargs, = (; obs_noise_cov = Γ * I(1)) +emulator = Emulator(mlt, iopairs; encoder_schedule = deepcopy(encoder_schedule), encoder_kwargs = encoder_kwargs) optimize_hyperparameters!(emulator) # (only needed for some Emulator packages) ``` diff --git a/docs/src/examples/emulators/lorenz_integrator_3d_3d.md b/docs/src/examples/emulators/lorenz_integrator_3d_3d.md index 87c40d14f..cea201361 100644 --- a/docs/src/examples/emulators/lorenz_integrator_3d_3d.md +++ b/docs/src/examples/emulators/lorenz_integrator_3d_3d.md @@ -90,10 +90,10 @@ end iopairs = PairedDataContainer(input, output) ``` We have several cases the user can play with, -```julia -cases = ["GP", "RF-scalar", "RF-scalar-diagin", "RF-svd-nonsep", "RF-nosvd-nonsep", "RF-nosvd-sep"] +```juliac +cases = ["GP", "RF-prior", "RF-lr-scalar", "rf-diag-scalar", "RF-lr-lr", "RF-nonsep"] ``` -Then, looping over the repeats, we first define some common hyperparamter optimization options for the `"RF-"` cases. In this case, the options are used primarily for diagnostics and acceleration (not required in general to solve this problem) +Then, looping over the repeats, we first define some common hyperparameter optimization options for the `"RF-"` cases. In this case, the options are used primarily for diagnostics and acceleration (not required in general to solve this problem) ```julia rf_optimizer_overrides = Dict( "verbose" => true, # output diagnostics from the hyperparameter optimizer @@ -113,7 +113,7 @@ mlt = GaussianProcess( noise_learn = false, # do not additionally learn a white kernel ) ``` -we also highlight the Vector Random Feature with nonseparable kernel (`RF-nosvd-nonsep`), this can natively handle multiple outputs without decorrelation of the output space. This kernel is a rank-3 representation with small nugget term. +we also highlight the Vector Random Feature with nonseparable kernel (`RF-nonsep`), this can natively handle multiple outputs without decorrelation of the output space. This kernel is a rank-3 representation with small nugget term. ```julia nugget = 1e-12 kernel_structure = NonseparableKernel(LowRankFactor(3, nugget)) @@ -129,9 +129,17 @@ mlt = VectorRandomFeatureInterface( ``` With machine learning tools specified, we build the emulator object ```julia -# decorrelate = true for `GP` -# decorrelate = false for `RF-nosvd-nonsep` -emulator = Emulator(mlt, iopairs; obs_noise_cov = Γy, decorrelate = decorrelate) +# e.g., for GP we decorrelate +encoder_schedule = (decorrelate_structure_mat(), "out") +encoder_kwargs = (; obs_noise_cov = Γy) + +emulator = Emulator( + mlt, + iopairs; + encoder_schedule = deepcopy(encoder_schedule), + encoder_kwargs = encoder_kwargs, +) + optimize_hyperparameters!(emulator) # some GP packages require this additional call ``` @@ -156,15 +164,15 @@ For one example fit ``` -### Random Feature Emulator (`RF-nosvd-nonsep`) +### Random Feature Emulator (`RF-nonsep`) For one example fit ```@raw html - - + + ``` and here are CDFs over 20 randomized trials of the random feature hyperparameter optimization ```@raw html - + ``` diff --git a/docs/src/examples/emulators/regression_2d_2d.md b/docs/src/examples/emulators/regression_2d_2d.md index 1dff8ed6c..2cd490681 100644 --- a/docs/src/examples/emulators/regression_2d_2d.md +++ b/docs/src/examples/emulators/regression_2d_2d.md @@ -34,21 +34,30 @@ To play with the hyperparameter optimization of RF, the optimizer options someti ```julia using CalibrateEmulateSample.EnsembleKalmanProcesses # Contains `DataMisfitController` ``` -We have 9 cases that the user can toggle or customize +We have several cases ofr (1) the machine learning tool, and (2) the type of data processing ```julia cases = [ "gp-skljl", "gp-gpjl", # Very slow prediction... - "rf-scalar", - "rf-svd-diag", - "rf-svd-nondiag", - "rf-nosvd-diag", - "rf-nosvd-nondiag", - "rf-svd-nonsep", - "rf-nosvd-nonsep", + "rf-lr-scalar", + "rf-lr-lr", + "rf-nonsep", ] + +encoder_names = [ + "no-proc", + "sample-struct-proc", # default processing + "sample-proc", + "struct-mat-proc", + "combined-proc", + "cca-proc", + "minmax", +] +# User chooses the mask vector +case_mask = [1, 3, 5] +encoder_mask = [1:7] ``` -The first two are for GP with either `ScikitLearn.jl` or `GaussianProcesses.jl` interface. The third is for the scalar RF interface, which most closely follows exactly replacing a GP. The rest are examples of vector RF with different types of data processing, (svd = same processing as scalar RF, nosvd = unprocessed) and different RF kernel structures in the output space of increasing complexity/flexibility (diag = Separable diagonal, nondiag = Separable nondiagonal, nonsep = nonseparable nondiagonal). +The first two are for GP with either `ScikitLearn.jl` or `GaussianProcesses.jl` interface. The third is for the scalar RF interface, which most closely follows exactly replacing a GP. The rest are examples of vector RF with different kernels. The different encodings use different matrices to do PCA, as well as a CCA option or a simple max-scaling option. We set up the learning problem specification, defining input and output dimensions, and number of data to train on, and the function `g` and the perturbed samples `y` with correlated additive noise ```julia @@ -89,8 +98,8 @@ We then build the emulators. An example for GP (`gp-skljl`) gppackage = SKLJL() # build a GP that learns an additional white noise kernel (along with the default RBF kernel) gaussian_process = GaussianProcess(gppackage, noise_learn = true) -# the data processing normalizes input data, and decorrelates output data with information from Σ -emulator = Emulator(gaussian_process, iopairs, obs_noise_cov = Σ, normalize_inputs = true) +encoder_kwargs = (; prior_cov = prior_cov, obs_noise_cov = Σ), +emulator = Emulator(gaussian_process, iopairs; encoder_schedule = deepcopy(encoder_schedule), encoder_kwargs = encoder_kwargs) ``` An example for scalar RF (`rf-scalar`) ```julia @@ -101,8 +110,7 @@ srfi = ScalarRandomFeatureInterface( kernel_structure = SeparableKernel(LowRankFactor(2, nugget), OneDimFactor()), optimizer_options = optimizer_options, ) -# the data processing normalizes input data, and decorrelates output data with information from Σ -emulator = Emulator(srfi, iopairs, obs_noise_cov = Σ, normalize_inputs = true) +emulator = Emulator(srfi, iopairs; encoder_schedule = deepcopy(encoder_schedule), encoder_kwargs = encoder_kwargs) ``` An example for vector RF (`rf-nosvd-nonsep`) ```julia @@ -114,8 +122,8 @@ vrfi = VectorRandomFeatureInterface( kernel_structure = NonseparableKernel(LowRankFactor(4, nugget)), optimizer_options = optimizer_options, ) -# the data processing normalizes input data, and does not decorrelate outputs -emulator = Emulator(vrfi, iopairs, obs_noise_cov = Σ, normalize_inputs = true, decorrelate = false) +encoder_kwargs = (; prior_cov = prior_cov, obs_noise_cov = Σ), +emulator = Emulator(vrfi, iopairs; encoder_schedule = deepcopy(encoder_schedule), encoder_kwargs = encoder_kwargs) ``` For RF and some GP packages, the training occurs during construction of the `Emulator`, however sometimes one must call an optimize step afterwards ``` diff --git a/docs/src/examples/lorenz_example.md b/docs/src/examples/lorenz_example.md index cdcd4a4f5..92c8987ad 100644 --- a/docs/src/examples/lorenz_example.md +++ b/docs/src/examples/lorenz_example.md @@ -139,20 +139,16 @@ The calibration stage must be run before the emulate-sample stages. The calibrat ## Emulate -Having run the `calibrate.jl` code to generate input-output pairs from parameters to data using `EnsembleKalmanProcesses.jl`, we will now run the Emulate and Sample stages (`emulate_sample.jl`). First, we need to define which machine learning model we will use for the emulation. We have 8 cases that the user can toggle or customize +Having run the `calibrate.jl` code to generate input-output pairs from parameters to data using `EnsembleKalmanProcesses.jl`, we will now run the Emulate and Sample stages (`emulate_sample.jl`). First, we need to define which machine learning model we will use for the emulation. We have a few cases that the user can toggle or customize ```julia cases = [ - "GP", # diagonalize, train scalar GP, assume diag inputs - "RF-scalar-diagin", # diagonalize, train scalar RF, assume diag inputs (most comparable to GP) - "RF-scalar", # diagonalize, train scalar RF, do not asume diag inputs - "RF-vector-svd-diag", - "RF-vector-svd-nondiag", - "RF-vector-nosvd-diag", - "RF-vector-nosvd-nondiag", - "RF-vector-svd-nonsep", + "gp-gpjl", # whiten and train scalar GP + "rf-lr-scalar", # whiten and train scalar-RF + "rf-lr-lr", # train vector-RF (separable kernel) + "rf-nonsep", # train vector-RF (nonseparable kernel) ] ``` -The first is for GP with `GaussianProcesses.jl` interface. The next two are for the scalar RF interface, which most closely follows exactly replacing a GP. The rest are examples of vector RF with different types of data processing, (svd = same processing as scalar RF, nosvd = unprocessed) and different RF kernel structures in the output space of increasing complexity/flexibility (diag = Separable diagonal, nondiag = Separable nondiagonal, nonsep = nonseparable nondiagonal). +The first is for GP with `GaussianProcesses.jl` interface. The next is for the scalar RF interface, which most closely follows exactly replacing a GP. The rest are examples of vector RF with different types of RF kernel structures in the output space of increasing complexity/flexibility. The example then loads the relevant training data that was constructed in the `calibrate.jl` call. ```julia @@ -184,10 +180,10 @@ mlt = GaussianProcess( ``` which calls `GaussianProcess.jl`. In this L96 example, since we focus on learning $F_s$ and $A$, we do not need to explicitly learn the noise, so `noise_learn = false`. -An example for scalar RF (`RF-scalar`) +An example for scalar RF (`rf-lr-scalar`) ```julia n_features = 100 -kernel_structure = SeparableKernel(LowRankFactor(2, nugget), OneDimFactor()) +kernel_structure = SeparableKernel(LowRankFactor(1, nugget), OneDimFactor()) mlt = ScalarRandomFeatureInterface( n_features, n_params, @@ -196,7 +192,7 @@ mlt = ScalarRandomFeatureInterface( optimizer_options = overrides, ) ``` -Optimizer options for `ScalarRandomFeature.jl` are provided throough `overrides` +Optimizer options for `ScalarRandomFeature.jl` are provided through `overrides` ```julia overrides = Dict( "verbose" => true, @@ -206,22 +202,26 @@ overrides = Dict( ) # we do not want termination, as our priors have relatively little interpretation ``` +We define the data processing. We scale the inputs with a robust scaling in each dimension. We the whiten the output space with `Γy`, and use a truncated PCA to retain 95% of the total variance): +```julia +retain_var = 0.95 +encoder_schedule = [(quartile_scale(), "in"), (decorrelate_structure_mat(retain_var = retain_var), "out")] +encoder_kwargs = (; obs_noise_cov = Γy) +``` We then build the emulator with the parameters as defined above ```julia +retain_var = 0.95 +encoder_schedule = [(quartile_scale(), "in"), (decorrelate_structure_mat(retain_var = retain_var), "out")] emulator = Emulator( mlt, input_output_pairs; - obs_noise_cov = Γy, - normalize_inputs = normalized, - standardize_outputs = standardize, - standardize_outputs_factors = norm_factor, - retained_svd_frac = retained_svd_frac, - decorrelate = decorrelate, + encoder_schedule=encoder_schedule, + encoder_kwargs=encoder_kwargs, ) ``` -For RF and some GP packages, the training occurs during construction of the Emulator, however sometimes one must call an optimize step afterwards +For RF and some GP packages, the training occurs during construction of the Emulator, however others require call an optimize step afterwards ```julia optimize_hyperparameters!(emulator) ``` diff --git a/docs/src/examples/sinusoid_example.md b/docs/src/examples/sinusoid_example.md index bfac67164..7b71e35a7 100644 --- a/docs/src/examples/sinusoid_example.md +++ b/docs/src/examples/sinusoid_example.md @@ -1,4 +1,4 @@ -# Sinusoid Example +# [Sinusoid Example](@id sinusoid-example) !!! info "How do I run this code?" The full code is found in the [`examples/`](https://github.com/CliMA/CalibrateEmulateSample.jl/tree/main/examples) directory of the github repository @@ -101,7 +101,7 @@ white_noise = MvNormal(zeros(dim_output), Γ) y_obs = [y1_true, y2_true] .+ rand(white_noise) println("Observations:", y_obs) ``` -This gives $y_{obs}=(6.15, 6.42)$. +This gives $y_{obs}=(6.01, 7.11)$. We can plot the true signal in black, the true observables in red and the noisy observables in blue. ![signal](../assets/sinusoid_true_vs_observed_signal.png) @@ -141,11 +141,11 @@ const PD = EKP.ParameterDistributions ``` We define prior distributions on the two parameters. For the amplitude, -we define a prior with mean 2 and standard deviation 1. It is +we define a prior with mean 4 and standard deviation 2. It is additionally constrained to be nonnegative. For the vertical shift we define -a Gaussian prior with mean 0 and standard deviation 5. +a Gaussian prior with mean 0 and standard deviation 8. ```julia -prior_u1 = PD.constrained_gaussian("amplitude", 2, 1, 0, Inf) +prior_u1 = PD.constrained_gaussian("amplitude", 4, 2, 0, Inf) prior_u2 = PD.constrained_gaussian("vert_shift", 0, 5, -Inf, Inf) prior = PD.combine_distributions([prior_u1, prior_u2]) # Plot priors @@ -153,7 +153,7 @@ p = plot(prior, fill = :lightgray) ``` ![prior](../assets/sinusoid_prior.png) -We now generate the initial ensemble and set up the EKI. +We now generate the initial ensemble and set up the EKI. We draw around 50 evaluations for illustration ```julia N_ensemble = 10 N_iterations = 5 @@ -181,14 +181,14 @@ Finally, we get the ensemble after the last iteration. This provides our estimat final_ensemble = EKP.get_ϕ_final(prior, ensemble_kalman_process) # Check that the ensemble mean is close to the theta_true -println("Ensemble mean: ", mean(final_ensemble, dims=2)) # [3.05, 6.37] +println("Ensemble mean: ", mean(final_ensemble, dims=2)) # [2.94, 7.08] println("True parameters: ", theta_true) # [3.0, 7.0] ``` | Parameter | Truth | EKI mean | | :---------------- | :------: | :----: | -| Amplitude | 3.0 | 3.05 | -| Vertical shift | 7.0 | 6.37 | +| Amplitude | 3.0 | 2.94 | +| Vertical shift | 7.0 | 7.08 | The EKI ensemble mean at the final iteration is close to the true parameters, which is good. We can also see how the ensembles evolve at each iteration in the plot below. @@ -248,8 +248,12 @@ See the [Gaussian process page](https://clima.github.io/CalibrateEmulateSample.j gppackage = Emulators.GPJL() gauss_proc = Emulators.GaussianProcess(gppackage, noise_learn = false) +# Create a data encoding to be performed on the input-output pairs +encoder_schedule = [(decorrelate_sample_cov(), "in"), (decorrelate_structure_mat(), "out")] +encoder_kwargs = (; obs_noise_cov = Γ,) + # Build emulator with data -emulator_gp = Emulator(gauss_proc, input_output_pairs, normalize_inputs = true, obs_noise_cov = Γ) +emulator_gp = Emulator(gauss_proc, input_output_pairs; encoder_schedule = encoder_schedule, encoder_kwargs = encoder_kwargs) optimize_hyperparameters!(emulator_gp) ``` For this simple example, we already know the observational noise `Γ=0.2*I`, so we set `noise_learn = false`. @@ -286,8 +290,9 @@ random_features = VectorRandomFeatureInterface( kernel_structure = kernel_structure, optimizer_options = optimizer_options, ) + emulator_random_features = - Emulator(random_features, input_output_pairs, normalize_inputs = true, obs_noise_cov = Γ, decorrelate = false) + Emulator(random_features, input_output_pairs; encoder_schedule = deepcopy(encoder_schedule), encoder_kwargs=deepcopy(encoder_kwargs)) optimize_hyperparameters!(emulator_random_features) ``` ### Emulator Validation @@ -362,7 +367,7 @@ We will provide the API with the observations, priors and our cheap emulator fro example we use the GP emulator. First, we need to find a suitable starting point, ideally one that is near the posterior distribution. We will use the final ensemble mean from EKI as this will increase the chance of acceptance near the start of the chain, and reduce burn-in time. ```julia init_sample = EKP.get_u_mean_final(ensemble_kalman_process) -println("initial parameters: ", init_sample) # (1.11, 6.37) +println("initial parameters: ", init_sample) # (1.07, 7.32) ``` Now, we can set up and carry out the MCMC starting from this point. @@ -379,10 +384,10 @@ chain = MarkovChainMonteCarlo.sample(mcmc, 100_000; stepsize = new_step, discard display(chain) ``` -| parameters | mean | std | -| :---------------- | :------: | :----: | -| amplitude | 1.1068 | 0.0943 | -| vert_shift | 6.3897 | 0.4601 | +| parameters | mean | std | +| :---------------- | :----: | :---: | +| amplitude | 1.10 | 0.073 | +| vert_shift | 7.07 | 0.46 | Note that these values are provided in the unconstrained space. The vertical shift seems reasonable, but the amplitude is not. This is because the amplitude is constrained to be @@ -394,7 +399,7 @@ constrained space and re-calculate these values. posterior = MarkovChainMonteCarlo.get_posterior(mcmc, chain) # Back to constrained coordinates constrained_posterior = Emulators.transform_unconstrained_to_constrained( - prior, MarkovChainMonteCarlo.get_distribution(posterior) + posterior, MarkovChainMonteCarlo.get_distribution(posterior) ) println("Amplitude mean: ", mean(constrained_posterior["amplitude"]), ", std: ", std(constrained_posterior["amplitude"])) println("Vertical shift mean: ", mean(constrained_posterior["vert_shift"]), ", std: ", std(constrained_posterior["vert_shift"])) @@ -402,10 +407,10 @@ println("Vertical shift mean: ", mean(constrained_posterior["vert_shift"]), ", s This gives: -| parameters | mean | std | -| :---------------- | :------: | :----: | -| amplitude | 3.0382 | 0.2880 | -| vert_shift | 6.3774 | 0.4586 | +| parameters | mean | std | +| :---------------- | :----: | :--: | +| amplitude | 3.01 | 0.21 | +| vert_shift | 7.07 | 0.46 | This is in agreement with the true $\theta=(3.0, 7.0)$ and with the observational covariance matrix we provided $\Gamma=0.2 * I$ (i.e., a standard deviation of approx. $0.45$). `CalibrateEmulateSample.jl` has built-in plotting recipes to help us visualize the prior and posterior distributions. Note that these are the @@ -433,10 +438,10 @@ $\theta_2$ below, with the marginal distributions on each axis. We can repeat the sampling method using the random features emulator instead of the Gaussian process and we find similar results: -| parameters | mean | std | -| :---------------- | :------: | :----: | -| amplitude | 3.3210 | 0.7216 | -| vert_shift | 6.3986 | 0.5098 | +| parameters | mean | std | +| :---------------- | :----: | :--: | +| amplitude | 2.99 | 0.21 | +| vert_shift | 7.06 | 4.89 | ![RF_2d_posterior](../assets/sinusoid_MCMC_hist_RF.png) diff --git a/docs/src/index.md b/docs/src/index.md index 28dc10da7..56bcc1e9e 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,10 +1,6 @@ # CalibrateEmulateSample.jl -`CalibrateEmulateSample.jl` solves parameter estimation problems using accelerated (and approximate) Bayesian inversion. - -The framework can be applied currently to learn: -- the joint distribution for a moderate numbers of parameters (<40), -- it is not inherently restricted to unimodal distributions. +`CalibrateEmulateSample.jl` solves parameter estimation problems using accelerated (and approximate) Bayesian inversion. The framework can be applied currently to learn the joint distribution for flexible numbers of parameters, it is not inherently restricted to unimodal distributions. Please see the [example below!](@ref inv-prob-front). It can be used with computer models that: - can be noisy or chaotic, @@ -19,14 +15,23 @@ y = \mathcal{G}(\theta) + \eta, ``` where the noise ``\eta`` is drawn from a $d$-dimensional Gaussian with distribution ``\mathcal{N}(0, \Gamma_y)``. +### Quick links + +- [How do I build prior distributions?](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/parameter_distributions/) +- [How do I build good observational noise covariances](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/observations/) +- [What ensemble size should I take? Which process should I use? What is the recommended configuration?](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/defaults/) +- [Where can I walk through the simple example?](@ref sinusoid-example) +- [What is the `EnsembleKalmanProcesses.jl` package?](@ref calibrate) +- [What are the recommendations/defaults for dimension reduction or data processing?](@ref data-proc) +- [How to I plot or interpret the posterior distribution?](@ref get-posterior) -### The inverse problem +### [The inverse problem](@id inv-prob-front) Given an observation ``y``, the computer model ``\mathcal{G}``, the observational noise ``\Gamma_y``, and some broad prior information on ``\theta``, we return the joint distribution of a data-informed distribution for "``\theta`` given ``y``". As the name suggests, `CalibrateEmulateSample.jl` breaks this problem into a sequence of three steps: calibration, emulation, and sampling. A comprehensive treatment of the calibrate-emulate-sample approach to Bayesian inverse problems can be found in [Cleary et al. (2020)](https://arxiv.org/pdf/2001.03689.pdf). -### The three steps of the algorithm: see our walkthrough of the [Sinusoid Example](@ref) +### The three steps of the algorithm: see our walkthrough of the [Sinusoid Example](@ref sinusoid-example) **Learn the vertical shift and amplitude of the signal given the noisy observation** ```@raw html diff --git a/docs/src/sample.md b/docs/src/sample.md index fdbc4d05f..2b686b590 100644 --- a/docs/src/sample.md +++ b/docs/src/sample.md @@ -8,30 +8,28 @@ The "sample" part of CES refers to exact sampling from the emulated posterior, i Carlo algorithm](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) (MCMC). Within this paradigm, we want to provide the flexibility to use multiple sampling algorithms; the approach we take is to use the general-purpose [AbstractMCMC.jl](https://turing.ml/dev/docs/for-developers/interface) API, provided by the [Turing.jl](https://turing.ml/dev/) probabilistic programming framework. -## User interface +## [User interface](@id sample-ui) We briefly outline an instance of how one sets up and uses MCMC within the CES package. The user first loads the MCMC module, and provides one of the Protocols (i.e. how one wishes to generate sampling proposals) - ```julia using CalibrateEmulateSample.MarkovChainMonteCarlo protocol = RWMHSampling() # Random-Walk algorithm -# protocol = pCNMHSampling() # preconditioned-Crank-Nicholson algorithm ``` -Then one builds the MCMC by providing the standard Bayesian ingredients (prior and data) from the calibrate stage, alongside the trained statistical emulator from the emulate stage: +More are listed below. One builds the MCMC by providing the standard Bayesian ingredients (prior and data) from the calibrate stage, alongside the trained statistical emulator from the emulate stage: ```julia mcmc = MCMCWrapper( protocol, - truth_sample, + observation_series, # or single observation prior, emulator; - init_params=mean_u_final, + init_params = mean_u_final, burnin=10_000, ) ``` -The keyword arguments `init_params` give a starting step of the chain (often taken to be the mean of the final iteration of calibrate stage), and a `burnin` gives a number of initial steps to be discarded when drawing statistics from the sampling method. +The keyword arguments `init_params` give a starting step of the chain (we recommend using the mean of the final iteration of calibrate stage), and a `burnin` gives a number of initial steps to be discarded when drawing statistics from the sampling method. -!!! note "for many samples" - If one has several samples of conditionally-independent data (that is, ``p({y_1,\dots,y_n}\mid\theta)`` is a product of ``\prod_i p(y_i\mid\theta)``), then one can feed in `truth_sample` as a vector of these samples, or a matrix with these samples as columns. The resulting sampler will evaluate the likelihood at all `y_i` for every sample step. +!!! note "Observations types" + If one has an `EnsembleKalmanProcesses` object (such as an `Observation` or `ObservationSeries`) simply pass this in as the second argument. If one has instead a sample of data simply pass this in, meanwhile a set of data samples can be passed in as a vector of vectors or matrix with these samples as columns. When multiple data is passed in (as an `ObservationSeries` or otherwise) the resulting sampler will assume the data is conditionally-independent (that is, ``p({y_1,\dots,y_n}\mid\theta)`` is a product of ``\prod_i p(y_i\mid\theta)``) and evaluate the likelihood at all `y_i` for every sample step. For good efficiency, one often needs to run MCMC with a problem-dependent step size. We provide a simple utility to help choose this. Here the optimizer runs short chains (of length `N`), and adjusts the step-size until the MCMC acceptance rate falls within an acceptable range, returning this step size. ```julia @@ -46,12 +44,20 @@ To generate ``10^5`` samples with a given step size (and optional random number chain = sample(rng, mcmc, 100_000; stepsize = new_step) display(chain) # gives diagnostics ``` -The return argument is stored in an `MCMCChains.Chains` object. To convert this back into a `ParameterDistribution` type (which contains e.g. the transformation maps) one can call +One can see our [examples](https://github.com/CliMA/CalibrateEmulateSample.jl/tree/main/examples) for other configurations! + +## [Interpretation and handling of the posterior samples](@id get-posterior) + +The posterior samples are under several wrappers when returned: +1. The first level (`MCMCChains.Chains`) is used to display e.g., MCMC diagnostics. +2. The second level (`Samples <: ParameterDistribution`) is used for easy plotting of histograms, and applying constraints. +3. The third level (`Dict` or `Array`) is the cleanest level of samples, and will be good for saving or applying other functions to without wrappers. + +To retrieve the `ParameterDistribution` object from the `MCMCChains.Chains` one can call ```julia posterior = get_posterior(mcmc, chain) -constrained_posterior = transform_unconstrained_to_constrained(prior, get_distribution(posterior)) +constrained_posterior = transform_unconstrained_to_constrained(posterior, get_distribution(posterior)) ``` - One can quickly plot the marginals of the prior and posterior distribution with ```julia using Plots @@ -63,6 +69,26 @@ or extract statistics of the (unconstrained) distribution with mean_posterior = mean(posterior) cov_posterior = cov(posterior) ``` +To convert the `ParameterDistribution` to `Dict` or `Array` types, in order to work with the parameters samples outside of any wrappers, call +```julia +# Dictionary-based: `Dict("param_name" => array_with_columns_as_samples)` +posterior_samples_dict = get_distribution(posterior) +constrained_posterior_samples_dict = transform_unconstrained_to_constrained(posterior, posterior_samples_dict) + +# Array-based (columns are samples) +posterior_samples_array = reduce(vcat, [posterior_samples_dict[n] for n in get_name(posterior)]) +constrained_posterior_samples_array = reduce(vcat, [constrained_posterior_samples_dict[n] for n in get_name(posterior)]) +``` + +## The Protocols + +For more details on the currently available protocols: +- `RWMHSampling()`: Random Walk Metropolis (standard). For example, see [Roberts, Gelman, Gilks (1997)](https://doi.org/10.1214/aoap/1034625254) +- `pCNMHSampling()`: preconditioned Crank Nicholsen (dimensionally-robust) For example, see [Cotter, Roberts, Stuart, White (2013)](https://www.jstor.org/stable/43288425) +- `BarkerSampling()` Barker proposal (derivative-based) For example, see [Livingstone, Zanella (2022)](https://doi.org/10.1111/rssb.12482). We currently have this working with two autodifferentiation packages, [`ForwardDiff.jl`](https://juliadiff.org/ForwardDiff.jl/stable/)(default) and [`ReverseDiff.jl`](https://juliadiff.org/ReverseDiff.jl/stable/). One can set the autodifferentiation protocols explicitly as `BarkerSampling{ReverseDiffProtocol}()`. `ReverseDiffProtocol` is typically slower for low dimensional problems, but scales better that `ForwardDiffProtocol` for higher dimensional settings. + +!!! warning "Differentiable samplers require differentiable emulators" + One caveat for differentiable samplers, is the lack of Julia emulator packages that support it for the `predict()` map, therefore one can only use e.g. `BarkerSampling` with [`AGPJL()`](@ref agpjl)-type emulators which currently only have limited support. # [Further details on the implementation](@id AbstractMCMC sampling API) diff --git a/examples/Lorenz/emulate_sample.jl b/examples/Lorenz/emulate_sample.jl index 1f13fc177..bbb8bf672 100644 --- a/examples/Lorenz/emulate_sample.jl +++ b/examples/Lorenz/emulate_sample.jl @@ -133,7 +133,6 @@ function main() ) end - # Standardize the output data # Get training points from the EKP iteration number in the second input term N_iter = min(max_iter, length(get_u(ekiobj)) - 1) # number of paired iterations taken from EKP diff --git a/examples/Sinusoid/calibrate.jl b/examples/Sinusoid/calibrate.jl index 330ef84be..740fe4797 100644 --- a/examples/Sinusoid/calibrate.jl +++ b/examples/Sinusoid/calibrate.jl @@ -68,8 +68,8 @@ data_save_directory = joinpath(example_directory, "output") # we define a prior with mean 2 and standard deviation 1. It is # additionally constrained to be nonnegative. For the vertical shift we define # a Gaussian prior with mean 0 and standard deviation 5. -prior_u1 = PD.constrained_gaussian("amplitude", 2, 1, 0, Inf) -prior_u2 = PD.constrained_gaussian("vert_shift", 0, 5, -Inf, Inf) +prior_u1 = PD.constrained_gaussian("amplitude", 4, 2, 0, Inf) +prior_u2 = PD.constrained_gaussian("vert_shift", 0, 8, -Inf, Inf) prior = PD.combine_distributions([prior_u1, prior_u2]) # Plot priors p = plot(prior, fill = :lightgray, rng = rng) diff --git a/examples/Sinusoid/emulate.jl b/examples/Sinusoid/emulate.jl index 06fdfca00..3cea02d9e 100644 --- a/examples/Sinusoid/emulate.jl +++ b/examples/Sinusoid/emulate.jl @@ -65,17 +65,25 @@ unconstrained_inputs = CES.Utilities.get_inputs(input_output_pairs) inputs = Emulators.transform_unconstrained_to_constrained(prior, unconstrained_inputs) outputs = CES.Utilities.get_outputs(input_output_pairs) +# Create a data encoding to be performed on the input-output pairs +encoder_schedule = [(decorrelate_sample_cov(), "in"), (decorrelate_structure_mat(), "out")] +encoder_kwargs = (; obs_noise_cov = Γ) + + # Gaussian process # We will set up a basic Gaussian process emulator using either ScikitLearn.jl or GaussianProcesses.jl. # See the Gaussian process page for more information and details about kernel choices. gppackage = Emulators.GPJL() gauss_proc = Emulators.GaussianProcess(gppackage, noise_learn = false) -# Create a data encoding to be performed on the input-output pairs -encoder_schedule = [(decorrelate_sample_cov(), "in"), (decorrelate_structure_mat(), "out")] # Build emulator with data -emulator_gp = Emulator(gauss_proc, input_output_pairs, (obs_cov_noise = Γ,); encoder_schedule = encoder_schedule) +emulator_gp = Emulator( + gauss_proc, + input_output_pairs; + encoder_schedule = deepcopy(encoder_schedule), + encoder_kwargs = deepcopy(encoder_kwargs), +) optimize_hyperparameters!(emulator_gp) # We have built the Gaussian process emulator and we can now use it for prediction. We will validate the emulator @@ -111,9 +119,13 @@ random_features = VectorRandomFeatureInterface( kernel_structure = kernel_structure, optimizer_options = optimizer_options, ) -encoder_schedule_rf = [(decorrelate_sample_cov(), "in"), (decorrelate_structure_mat(), "out")] -emulator_random_features = - Emulator(random_features, input_output_pairs, (obs_cov_noise = Γ,); encoder_schedule = encoder_schedule_rf) + +emulator_random_features = Emulator( + random_features, + input_output_pairs; + encoder_schedule = deepcopy(encoder_schedule), + encoder_kwargs = deepcopy(encoder_kwargs), +) optimize_hyperparameters!(emulator_random_features) @@ -123,8 +135,8 @@ optimize_hyperparameters!(emulator_random_features) # In more complex systems, we would have limited data to validate emulator performance with. # Set up mesh to cover parameter space N_grid = 50 -amp_range = range(0.6, 4, length = N_grid) -vshift_range = range(-6, 10, length = N_grid) +amp_range = range(0.5, 6, length = N_grid) +vshift_range = range(-15, 15, length = N_grid) function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T} m, n = length(vy), length(vx) @@ -160,8 +172,8 @@ rf_std_grid = reshape(permutedims(reduce(vcat, [x' for x in rf_std]), (2, 1)), ( ## Plot # First, we will plot the ground truth. We have 2 parameters and 2 outputs, so we will create a contour plot # for each output to show how they vary against the two inputs. -range_clims = (0, 8) -mean_clims = (-6, 10) +range_clims = (0, 12) +mean_clims = (-15, 15) p1 = contour( amp_range, @@ -271,7 +283,7 @@ savefig(p, joinpath(data_save_directory, "sinusoid_RF_emulator_contours.png")) # Next, we plot uncertainty estimates from the GP and RF emulators # Plot GP std estimates -range_std_clims = (0, 2) +range_std_clims = (0, 1) mean_std_clims = (0, 1) p1 = contour( @@ -341,8 +353,8 @@ savefig(p, joinpath(data_save_directory, "sinusoid_RF_emulator_std_contours.png" # Finally, we should validate how accurate the emulators are by looking at the absolute difference between emulator # predictions and the ground truth. gp_diff_grid = abs.(gp_grid - g_true_grid) -range_diff_clims = (0, 1) -mean_diff_clims = (0, 1) +range_diff_clims = (0, 2) +mean_diff_clims = (0, 2) p1 = contour( amp_range, vshift_range, diff --git a/src/Utilities/decorrelator.jl b/src/Utilities/decorrelator.jl index aa94266f2..6c5ec24c6 100644 --- a/src/Utilities/decorrelator.jl +++ b/src/Utilities/decorrelator.jl @@ -56,7 +56,7 @@ struct Decorrelator{VV1, VV2, VV3, FT, NT <: NamedTuple, AS <: AbstractString} < n_totvar_samples::Int "maximum dimension of subspace for `retain_var < 1`. The search may become expensive at large ranks, and therefore can be cut-off in this way" max_rank::Int - "when `retain_var = 1`, the `psvd` algorithm from `LowRankApprox.jl` is used to decorrelate the space. here, kwargs can be passed in as a NamedTuple" + "when `retain_var = 1`, the `psvd` algorithm from `LowRankApprox.jl` is used to decorrelate the space. here, kwargs can be passed in as a `NamedTuple`" psvd_kwargs::NT "Switch to choose what form of matrix to use to decorrelate the data" decorrelate_with::AS @@ -73,6 +73,9 @@ Constructs the `Decorrelator` struct. Users can add optional keyword arguments: - `"structure_mat"`, see [`decorrelate_structure_mat`](@ref) - `"sample_cov"`, see [`decorrelate_sample_cov`](@ref) - `"combined"`, sums the `"sample_cov"` and `"structure_mat"` matrices +- `n_totvar_samples`[=`500`]: when retain_var < 1, number of samples to estimate the total variance for performing truncation. +- `max_rank`[=`100`]: for `retain_var < 1`, the maximum dimension of subspace when using an the `tsvd` algorithm from `TSVD.jl`. +- `psvd_kwargs` [= `(; rtol = 1e-3)`]: for `retain_var = 1`, the `psvd` algorithm from `LowRankApprox.jl` is used to decorrelate the space. kwargs can be passed in as a `NamedTuple` """ decorrelate(; retain_var::FT = Float64(1.0), @@ -98,6 +101,9 @@ $(TYPEDSIGNATURES) Constructs the `Decorrelator` struct, setting decorrelate_with = "sample_cov". Encoding data with this will ensure that the distribution of data samples after encoding will be `Normal(0,I)`. One can additionally add keywords: - `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained +- `n_totvar_samples`[=`500`]: when retain_var < 1, number of samples to estimate the total variance for performing truncation. +- `max_rank`[=`100`]: for `retain_var < 1`, the maximum dimension of subspace when using an the `tsvd` algorithm from `TSVD.jl`. +- `psvd_kwargs` [= `(; rtol = 1e-3)`]: for `retain_var = 1`, the `psvd` algorithm from `LowRankApprox.jl` is used to decorrelate the space. kwargs can be passed in as a `NamedTuple` """ decorrelate_sample_cov(; retain_var::FT = Float64(1.0), @@ -121,6 +127,9 @@ $(TYPEDSIGNATURES) Constructs the `Decorrelator` struct, setting decorrelate_with = "structure_mat". This encoding will transform a provided structure matrix into `I`. One can additionally add keywords: - `retain_var`[=`1.0`]: to project onto the leading singular vectors such that `retain_var` variance is retained +- `n_totvar_samples`[=`500`]: when retain_var < 1, number of samples to estimate the total variance for performing truncation. +- `max_rank`[=`100`]: for `retain_var < 1`, the maximum dimension of subspace when using an the `tsvd` algorithm from `TSVD.jl`. +- `psvd_kwargs` [= `(; rtol = 1e-3)`]: for `retain_var = 1`, the `psvd` algorithm from `LowRankApprox.jl` is used to decorrelate the space. kwargs can be passed in as a `NamedTuple` """ decorrelate_structure_mat(; retain_var::FT = Float64(1.0),