CliMA · odunbar · Jan 10, 2026 · Dec 19, 2025 · Dec 19, 2025 · Dec 19, 2025
diff --git 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
@@ -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
diff --git a/docs/src/assets/GFunction_sens_RF-scalar_3.png b/docs/src/assets/GFunction_sens_RF-scalar_3.png
diff --git a/docs/src/assets/GFunction_slices_RF-scalar_10.png b/docs/src/assets/GFunction_slices_RF-scalar_10.png
diff --git a/docs/src/assets/GFunction_slices_RF-scalar_3.png b/docs/src/assets/GFunction_slices_RF-scalar_3.png
diff --git a/docs/src/assets/GP_l63_attr.png b/docs/src/assets/GP_l63_attr.png
diff --git a/docs/src/assets/GP_l63_cdfs.png b/docs/src/assets/GP_l63_cdfs.png
diff --git a/docs/src/assets/GP_l63_pdf.png b/docs/src/assets/GP_l63_pdf.png
diff --git a/docs/src/assets/GP_l63_test.png b/docs/src/assets/GP_l63_test.png
diff --git a/docs/src/assets/RF-nonsep_l63_attr.png b/docs/src/assets/RF-nonsep_l63_attr.png
diff --git a/docs/src/assets/RF-nonsep_l63_cdfs.png b/docs/src/assets/RF-nonsep_l63_cdfs.png
diff --git a/docs/src/assets/RF-nonsep_l63_pdf.png b/docs/src/assets/RF-nonsep_l63_pdf.png
diff --git a/docs/src/assets/RF-nonsep_l63_test.png b/docs/src/assets/RF-nonsep_l63_test.png
diff --git a/docs/src/assets/RF-nosvd-nonsep_l63_attr.png b/docs/src/assets/RF-nosvd-nonsep_l63_attr.png
diff --git a/docs/src/assets/RF-nosvd-nonsep_l63_cdfs.png b/docs/src/assets/RF-nosvd-nonsep_l63_cdfs.png
diff --git a/docs/src/assets/RF-nosvd-nonsep_l63_pdf.png b/docs/src/assets/RF-nosvd-nonsep_l63_pdf.png
diff --git a/docs/src/assets/RF-nosvd-nonsep_l63_test.png b/docs/src/assets/RF-nosvd-nonsep_l63_test.png
diff --git a/docs/src/assets/sinusoid_GP_emulator_contours.png b/docs/src/assets/sinusoid_GP_emulator_contours.png
diff --git a/docs/src/assets/sinusoid_GP_errors_contours.png b/docs/src/assets/sinusoid_GP_errors_contours.png
diff --git a/docs/src/assets/sinusoid_MCMC_hist_GP.png b/docs/src/assets/sinusoid_MCMC_hist_GP.png
diff --git a/docs/src/assets/sinusoid_MCMC_hist_RF.png b/docs/src/assets/sinusoid_MCMC_hist_RF.png
diff --git a/docs/src/assets/sinusoid_RF_emulator_contours.png b/docs/src/assets/sinusoid_RF_emulator_contours.png
diff --git a/docs/src/assets/sinusoid_RF_errors_contours.png b/docs/src/assets/sinusoid_RF_errors_contours.png
diff --git a/docs/src/assets/sinusoid_eki_pairs.png b/docs/src/assets/sinusoid_eki_pairs.png
diff --git a/docs/src/assets/sinusoid_groundtruth_contours.png b/docs/src/assets/sinusoid_groundtruth_contours.png
diff --git a/docs/src/assets/sinusoid_posterior_GP.png b/docs/src/assets/sinusoid_posterior_GP.png
diff --git a/docs/src/assets/sinusoid_posterior_RF.png b/docs/src/assets/sinusoid_posterior_RF.png
diff --git a/docs/src/assets/sinusoid_prior.png b/docs/src/assets/sinusoid_prior.png
diff --git a/docs/src/assets/sinusoid_true_vs_observed_signal.png b/docs/src/assets/sinusoid_true_vs_observed_signal.png
diff --git 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
@@ -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
@@ -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
@@ -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)
 ```