diff --git a/docs/src/data_processing.md b/docs/src/data_processing.md
index e1db4e328..662bcad95 100644
--- a/docs/src/data_processing.md
+++ b/docs/src/data_processing.md
@@ -7,8 +7,9 @@ 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. 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.
+- Users can predict to/from encoded space with `predict(...; encode=X)` with `X` being `"in"`, `"out"`, or `"in_and_out"` if needed. 
 
-An external API is also available using the `encode_data`, and `encode_structure_matrix` methods if needed.
+An external API is also available using the `encode_data`, and `encode_structure_matrix` methods if needed to encode new data, or new covariance-type matrices.
 
 ## Defaults and recommendations
 ### Default schedule (no explicit dimension reduction)
diff --git a/docs/src/emulate.md b/docs/src/emulate.md
index f53cb4f42..bfb2d7f42 100644
--- a/docs/src/emulate.md
+++ b/docs/src/emulate.md
@@ -33,9 +33,16 @@ optimize_hyperparameters!(emulator)
 For some machine learning packages however, this may be completed during construction automatically, and for others this will not. If automatic construction took place, the `optimize_hyperparameters!` line does not perform any new task, so may be safely called. In the Lorenz example, this line learns the hyperparameters of the Gaussian process, which depend on the choice of [kernel](https://clima.github.io/CalibrateEmulateSample.jl/dev/GaussianProcessEmulator/#kernels), and the choice of GP package.
 Predictions at new inputs can then be made using
 ```julia
-y, cov = Emulator.predict(emulator, new_inputs)
+em_mean, em_cov = Emulator.predict(emulator, new_inputs)
+```
+This returns both a mean value and a covariance. The emulator is subject to encoding (see [Data processing and Dimension reduction](@ref data-proc)), and so we provide the `encode` and `add_obs_noise_cov` to enable users to predict in different spaces, and with different inflation.
+```julia
+# produce output in encoded space
+em_mean_enc, em_cov_enc = Emulator.predict(emulator, new_inputs; encode="out")
+
+# given encoded inputs, produce outputs in real space, and inflate the emulator uncertainty with observational noise
+em_mean, em_and_obs_noise_cov = Emulator.predict(emulator, new_encoded_inputs; encode="in", add_obs_noise_cov=true) 
 ```
-This returns both a mean value and a covariance.
 
 ## [Modular interface](@id modular-interface)
 
diff --git a/docs/src/examples/Cloudy_example.md b/docs/src/examples/Cloudy_example.md
index 0c227f57d..e043974f7 100755
--- a/docs/src/examples/Cloudy_example.md
+++ b/docs/src/examples/Cloudy_example.md
@@ -417,8 +417,7 @@ on them).
 # Check how well the emulator predicts on the true parameters
 y_mean, y_var = Emulators.predict(
     emulator,
-    reshape(θ_true, :, 1);
-    transform_to_real = true
+    reshape(θ_true, :, 1)
 )
 
 println("Emulator ($(case)) prediction on true parameters: ")
diff --git a/docs/src/examples/emulators/global_sens_analysis.md b/docs/src/examples/emulators/global_sens_analysis.md
index 3bb794c16..210e41c2d 100644
--- a/docs/src/examples/emulators/global_sens_analysis.md
+++ b/docs/src/examples/emulators/global_sens_analysis.md
@@ -127,11 +127,11 @@ optimize_hyperparameters!(emulator) # (only needed for some Emulator packages)
 
 ### Results and plots
 
-We validate the emulator by evaluating it on the entire Sobol sequence, and calculating the Sobol indices (presenting mean and std if using multiple repeats.
+We validate the emulator by evaluating it on the entire Sobol sequence, and calculating the Sobol indices (presenting mean and std if using multiple repeats).
 
 ```julia
 # predict on all Sobol points with emulator (example)    
-y_pred, y_var = predict(emulator, samples', transform_to_real = true)
+y_pred, y_var = predict(emulator, samples')
 
 # obtain emulated Sobol indices
 result_pred = analyze(data, y_pred')
diff --git a/docs/src/examples/emulators/lorenz_integrator_3d_3d.md b/docs/src/examples/emulators/lorenz_integrator_3d_3d.md
index cea201361..249e24a7e 100644
--- a/docs/src/examples/emulators/lorenz_integrator_3d_3d.md
+++ b/docs/src/examples/emulators/lorenz_integrator_3d_3d.md
@@ -151,7 +151,7 @@ u_test_tmp = zeros(3, length(xspan_test))
 u_test_tmp[:, 1] = sol_test.u[1] # initialize at the final-time solution of the training period
 
 for i in 1:(length(xspan_test) - 1)
-    rf_mean, _ = predict(emulator, u_test_tmp[:, i:i], transform_to_real = true) # 3x1 matrix
+    rf_mean, _ = predict(emulator, u_test_tmp[:, i:i]) # 3x1 matrix
     u_test_tmp[:, i + 1] = rf_mean
 end
 ```
diff --git a/docs/src/examples/emulators/regression_2d_2d.md b/docs/src/examples/emulators/regression_2d_2d.md
index 2cd490681..6b0af4d2a 100644
--- a/docs/src/examples/emulators/regression_2d_2d.md
+++ b/docs/src/examples/emulators/regression_2d_2d.md
@@ -139,9 +139,9 @@ x2 = range(0.0, stop = 2 * π, length = n_pts)
 X1, X2 = meshgrid(x1, x2)
 inputs = permutedims(hcat(X1[:], X2[:]), (2, 1))
 ```
-We predict using the emulators at the new inputs, and `transform_to_real` inverts the data processing back to physical values
+We predict using the emulators at the new inputs, 
 ```julia
-em_mean, em_cov = predict(emulator, inputs, transform_to_real = true)
+em_mean, em_cov = predict(emulator, inputs)
 ```
 We then plot the predicted mean and pointwise variances, and calculate the errors from the three highlighted cases:
 
diff --git a/docs/src/examples/lorenz_example.md b/docs/src/examples/lorenz_example.md
index 92c8987ad..6939b3684 100644
--- a/docs/src/examples/lorenz_example.md
+++ b/docs/src/examples/lorenz_example.md
@@ -230,8 +230,8 @@ The emulator is checked for accuracy by evaluating its predictions on the true p
 ```julia
 # Check how well the Gaussian Process regression predicts on the
 # true parameters
-y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1), transform_to_real = true)
-y_mean_test, y_var_test = Emulators.predict(emulator, get_inputs(input_output_pairs_test), transform_to_real = true)
+y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1))
+y_mean_test, y_var_test = Emulators.predict(emulator, get_inputs(input_output_pairs_test))
 
 println("ML prediction on true parameters: ")
 println(vec(y_mean))
diff --git a/examples/Cloudy/Cloudy_emulate_sample.jl b/examples/Cloudy/Cloudy_emulate_sample.jl
index 4c1a2b514..b71b5f687 100644
--- a/examples/Cloudy/Cloudy_emulate_sample.jl
+++ b/examples/Cloudy/Cloudy_emulate_sample.jl
@@ -181,8 +181,8 @@ function main()
         optimize_hyperparameters!(emulator)
 
         # Check how well the emulator predicts on the true parameters
-        y_mean, y_var = Emulators.predict(emulator, reshape(θ_true, :, 1); transform_to_real = true)
-        y_mean_test, y_var_test = Emulators.predict(emulator, get_inputs(test_pairs); transform_to_real = true)
+        y_mean, y_var = Emulators.predict(emulator, reshape(θ_true, :, 1))
+        y_mean_test, y_var_test = Emulators.predict(emulator, get_inputs(test_pairs))
         println("Emulator ($(case)) prediction on true parameters: ")
         println(vec(y_mean))
         println("true data: ")
diff --git a/examples/Darcy/emulate_sample.jl b/examples/Darcy/emulate_sample.jl
index b7874a3c2..6bad3788f 100644
--- a/examples/Darcy/emulate_sample.jl
+++ b/examples/Darcy/emulate_sample.jl
@@ -105,9 +105,8 @@ function main()
         # Check how well the Gaussian Process regression predicts on the
         # true parameters
         #if retained_svd_frac==1.0
-        y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1), transform_to_real = true)
-        y_mean_test, y_var_test =
-            Emulators.predict(emulator, get_inputs(input_output_pairs_test), transform_to_real = true)
+        y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1))
+        y_mean_test, y_var_test = Emulators.predict(emulator, get_inputs(input_output_pairs_test))
 
         println("ML prediction on true parameters: ")
         println(vec(y_mean))
diff --git a/examples/EDMF_data/emulator-rank-test.jl b/examples/EDMF_data/emulator-rank-test.jl
index 9817de632..82ab02f7f 100644
--- a/examples/EDMF_data/emulator-rank-test.jl
+++ b/examples/EDMF_data/emulator-rank-test.jl
@@ -268,7 +268,7 @@ function main()
             # error of emulator on the training points: (no denoised training data)
             train_err_tmp = [0.0]
             for i in 1:size(train_inputs, 2)
-                train_mean, _ = Emulators.predict(emulator, train_inputs[:, i:i], transform_to_real = true) # 3x1
+                train_mean, _ = Emulators.predict(emulator, train_inputs[:, i:i]) # 3x1
                 train_err_tmp[1] += norm(train_mean - train_outputs[:, i])
             end
             train_err[rank_id, rep_idx] = 1 / size(train_inputs, 2) * train_err_tmp[1]
@@ -276,7 +276,7 @@ function main()
             # error of emulator on test points: 
             test_err_tmp = [0.0]
             for i in 1:size(test_inputs, 2)
-                test_mean, _ = Emulators.predict(emulator, test_inputs[:, i:i], transform_to_real = true) # 3x1
+                test_mean, _ = Emulators.predict(emulator, test_inputs[:, i:i]) # 3x1
                 test_err_tmp[1] += norm(test_mean - test_outputs[:, i])
             end
             test_err[rank_id, rep_idx] = 1 / size(test_inputs, 2) * test_err_tmp[1]
diff --git a/examples/Emulator/G-function/emulate-test-n-features.jl b/examples/Emulator/G-function/emulate-test-n-features.jl
index ecb324a94..5591038b8 100644
--- a/examples/Emulator/G-function/emulate-test-n-features.jl
+++ b/examples/Emulator/G-function/emulate-test-n-features.jl
@@ -162,11 +162,11 @@ function main()
 
             # errors:
             # training error
-            y_pred, y_var = predict(emulator, get_inputs(iopairs), transform_to_real = true)
+            y_pred, y_var = predict(emulator, get_inputs(iopairs))
             train_err[f_idx, rep_idx] = sqrt(sum((y_pred - get_outputs(iopairs)) .^ 2)) / n_train_pts
 
             # predict on all test points with emulator (example)    
-            y_pred, y_var = predict(emulator, samples', transform_to_real = true) #predict on all points
+            y_pred, y_var = predict(emulator, samples') #predict on all points
             ind_test = ind_total[(n_train_pts + 1):end]
             test_err[f_idx, rep_idx] = sqrt(sum((y_pred[ind_test] - y[ind_test]) .^ 2)) / length(ind_test)
 
diff --git a/examples/Emulator/G-function/emulate.jl b/examples/Emulator/G-function/emulate.jl
index 6c4afd02e..d3d8820b4 100644
--- a/examples/Emulator/G-function/emulate.jl
+++ b/examples/Emulator/G-function/emulate.jl
@@ -159,7 +159,7 @@ function main()
 
         @info "statistics of training time for case $(case): \n mean(s): $(mean(times[1:rep_idx])) \n var(s) : $(var(times[1:rep_idx]))"
         # predict on all Sobol points with emulator (example)    
-        y_pred, y_var = predict(emulator, samples', transform_to_real = true)
+        y_pred, y_var = predict(emulator, samples')
 
         # obtain emulated Sobol indices
         result_pred = analyze(data, y_pred')
diff --git a/examples/Emulator/Ishigami/emulate.jl b/examples/Emulator/Ishigami/emulate.jl
index 6136590ac..0be8b7922 100644
--- a/examples/Emulator/Ishigami/emulate.jl
+++ b/examples/Emulator/Ishigami/emulate.jl
@@ -138,7 +138,7 @@ function main()
             push!(opt_diagnostics, diag_tmp)
         end
         # predict on all Sobol points with emulator (example)    
-        y_pred, y_var = predict(emulator, samples', transform_to_real = true)
+        y_pred, y_var = predict(emulator, samples')
 
         # obtain emulated Sobol indices
         result_pred = analyze(data, y_pred')
diff --git a/examples/Emulator/L63/emulate.jl b/examples/Emulator/L63/emulate.jl
index 026e3c51a..49a57f176 100644
--- a/examples/Emulator/L63/emulate.jl
+++ b/examples/Emulator/L63/emulate.jl
@@ -214,13 +214,13 @@ function main()
         u_test_tmp[:, 1] = sol_test.u[1]
 
         for i in 1:(length(xspan_test) - 1)
-            rf_mean, _ = predict(emulator, u_test_tmp[:, i:i], transform_to_real = true) # 3x1 matrix
+            rf_mean, _ = predict(emulator, u_test_tmp[:, i:i]) # 3x1 matrix
             u_test_tmp[:, i + 1] = rf_mean
         end
 
         train_err_tmp = [0.0]
         for i in 1:size(input, 2)
-            train_mean, _ = predict(emulator, input[:, i:i], transform_to_real = true) # 3x1
+            train_mean, _ = predict(emulator, input[:, i:i]) # 3x1
             train_err_tmp[1] += norm(train_mean - output[:, i])
         end
         println("normalized L^2 error on training data:", 1 / size(input, 2) * train_err_tmp[1])
@@ -229,7 +229,7 @@ function main()
         u_hist_tmp[:, 1] = sol_hist.u[1]  # start at end of previous sim
 
         for i in 1:(length(xspan_hist) - 1)
-            rf_mean, _ = predict(emulator, u_hist_tmp[:, i:i], transform_to_real = true) # 3x1 matrix
+            rf_mean, _ = predict(emulator, u_hist_tmp[:, i:i]) # 3x1 matrix
             u_hist_tmp[:, i + 1] = rf_mean
         end
 
diff --git a/examples/Emulator/L63/emulate_diff-rank-test.jl b/examples/Emulator/L63/emulate_diff-rank-test.jl
index 38c896a1f..3d8c54fe0 100644
--- a/examples/Emulator/L63/emulate_diff-rank-test.jl
+++ b/examples/Emulator/L63/emulate_diff-rank-test.jl
@@ -211,14 +211,14 @@ function main()
 
             # predict sequentially i -> i+1
             for i in 1:(length(xspan_test) - 1)
-                rf_mean, _ = predict(emulator, u_test_tmp[:, i:i], transform_to_real = true) # 3x1 matrix
+                rf_mean, _ = predict(emulator, u_test_tmp[:, i:i]) # 3x1 matrix
                 u_test_tmp[:, i + 1] = rf_mean
             end
 
             # training error i -> o
             train_err_tmp = [0.0]
             for i in 1:size(input, 2)
-                train_mean, _ = predict(emulator, input[:, i:i], transform_to_real = true) # 3x1
+                train_mean, _ = predict(emulator, input[:, i:i]) # 3x1
                 train_err_tmp[1] += norm(train_mean - output[:, i])
             end
             train_err[rank_id, rep_idx] = 1 / size(input, 2) * train_err_tmp[1]
@@ -226,7 +226,7 @@ function main()
             # test error i -> o
             test_err_tmp = [0.0]
             for i in 1:(length(xspan_test) - 1)
-                test_mean, _ = predict(emulator, reshape(sol_test.u[i], :, 1), transform_to_real = true) # 3x1 matrix
+                test_mean, _ = predict(emulator, reshape(sol_test.u[i], :, 1)) # 3x1 matrix
                 test_err_tmp[1] += norm(test_mean[:] - sol_test.u[i + 1])
             end
             test_err[rank_id, rep_idx] = 1 / (length(xspan_test) - 1) * test_err_tmp[1]
@@ -237,7 +237,7 @@ function main()
             u_hist_tmp[:, 1] = sol_hist.u[1]  # start at end of previous sim
 
             for i in 1:(length(xspan_hist) - 1)
-                rf_mean, _ = predict(emulator, u_hist_tmp[:, i:i], transform_to_real = true) # 3x1 matrix
+                rf_mean, _ = predict(emulator, u_hist_tmp[:, i:i]) # 3x1 matrix
                 u_hist_tmp[:, i + 1] = rf_mean
             end
 
diff --git a/examples/Emulator/Regression_2d_2d/compare_regression.jl b/examples/Emulator/Regression_2d_2d/compare_regression.jl
index 2becbeeea..e7be1701d 100644
--- a/examples/Emulator/Regression_2d_2d/compare_regression.jl
+++ b/examples/Emulator/Regression_2d_2d/compare_regression.jl
@@ -206,7 +206,7 @@ function main()
             # Input for predict has to be of size  input_dim x N_samples
             inputs = permutedims(hcat(X1[:], X2[:]), (2, 1))
 
-            em_mean, em_cov = predict(emulator, inputs, transform_to_real = true)
+            em_mean, em_cov = predict(emulator, inputs)
             println("end predictions at ", n_pts * n_pts, " points")
 
             g1_true = g1(inputs)
diff --git a/examples/GCM/emulate_sample_script.jl b/examples/GCM/emulate_sample_script.jl
index 370996d15..f861b2d8d 100644
--- a/examples/GCM/emulate_sample_script.jl
+++ b/examples/GCM/emulate_sample_script.jl
@@ -197,7 +197,7 @@ function main()
             new_input = [log(0.31735951644387783 / (1 - 0.31735951644387783)) log(90632.50269636544)]' # random parameter value ("far" from truth
         end
 
-        pred_mean, pred_cov = predict(emulator, new_input, transform_to_real = true)
+        pred_mean, pred_cov = predict(emulator, new_input)
         pred_sd = sqrt.([max(10 * eps(), pred_cov[1][i, i]) for i in 1:size(pred_cov[1], 1)])
 
 
@@ -223,7 +223,7 @@ function main()
     end
 
     plot_input = [log(0.7 / 0.3) log(7200)]' # physical parameter value (at truth)
-    plot_mean, plot_cov = predict(emulator, plot_input, transform_to_real = true)
+    plot_mean, plot_cov = predict(emulator, plot_input)
     plot_sd = sqrt.([max(10 * eps(), plot_cov[1][i, i]) for i in 1:size(plot_cov[1], 1)])
 
 
diff --git a/examples/Lorenz/emulate_sample.jl b/examples/Lorenz/emulate_sample.jl
index bbb8bf672..b7466463c 100644
--- a/examples/Lorenz/emulate_sample.jl
+++ b/examples/Lorenz/emulate_sample.jl
@@ -181,9 +181,8 @@ function main()
         # Check how well the Gaussian Process regression predicts on the
         # true parameters
         #if retained_svd_frac==1.0
-        y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1), transform_to_real = true)
-        y_mean_test, y_var_test =
-            Emulators.predict(emulator, get_inputs(input_output_pairs_test), transform_to_real = true)
+        y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1))
+        y_mean_test, y_var_test = Emulators.predict(emulator, get_inputs(input_output_pairs_test))
 
         println("ML prediction on true parameters: ")
         println(vec(y_mean))
diff --git a/examples/Lorenz/emulate_sample_spatial_dep.jl b/examples/Lorenz/emulate_sample_spatial_dep.jl
index 9cffec89e..aa9aa758e 100644
--- a/examples/Lorenz/emulate_sample_spatial_dep.jl
+++ b/examples/Lorenz/emulate_sample_spatial_dep.jl
@@ -177,9 +177,8 @@ function main()
 
         # Check how well the Gaussian Process regression predicts on the
         # true parameters
-        y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1), transform_to_real = true)
-        y_mean_test, y_var_test =
-            Emulators.predict(emulator, get_inputs(input_output_pairs_test), transform_to_real = true)
+        y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1))
+        y_mean_test, y_var_test = Emulators.predict(emulator, get_inputs(input_output_pairs_test))
 
         println("ML prediction on true parameters: ")
         println(vec(y_mean))
diff --git a/examples/Sinusoid/emulate.jl b/examples/Sinusoid/emulate.jl
index 74d5bbe29..9b40e104e 100644
--- a/examples/Sinusoid/emulate.jl
+++ b/examples/Sinusoid/emulate.jl
@@ -168,9 +168,9 @@ g_true_grid = reshape(g_true, (output_dim, N_grid, N_grid))
 # Next, predict with emulators. We first need to transform to the unconstrained space.
 input_grid_unconstrained = Emulators.transform_constrained_to_unconstrained(prior, input_grid)
 
-gp_mean, gp_cov = Emulators.predict(emulator_gp, input_grid_unconstrained, transform_to_real = true)
-rf_mean, rf_cov = Emulators.predict(emulator_rf, input_grid_unconstrained, transform_to_real = true)
-fm_mean, fm_cov = Emulators.predict(emulator_fm, input_grid_unconstrained; transform_to_real = true)
+gp_mean, gp_cov = Emulators.predict(emulator_gp, input_grid_unconstrained)
+rf_mean, rf_cov = Emulators.predict(emulator_rf, input_grid_unconstrained)
+fm_mean, fm_cov = Emulators.predict(emulator_fm, input_grid_unconstrained)
 
 
 # Reshape into (output_dim x 50 x 50) grid
diff --git a/src/Emulator.jl b/src/Emulator.jl
index 1caf82c89..a4507f078 100644
--- a/src/Emulator.jl
+++ b/src/Emulator.jl
@@ -356,10 +356,10 @@ Makes a prediction using the emulator on new inputs (each new inputs given as da
 
 Keyword args
 - `encode` [=`nothing`]: For the input encoder `Eᵢ`, and output decoder `Dₒ` stored in the emulator, we have learnt a predict method method `G` in the encoded space. Interpret the keyword as follows:
-    - nothing     : applies Dₒ∘G∘Eᵢ(x) (nothing is encoded) - most common for user interaction
-    - "in"        : applies Dₒ∘G(z) (the inputs are provided as encoded (z=Eᵢx))
-    - "out"       : applies G∘Eᵢ(x) (the outputs are returned as encoded)
-    - "in_and_out": applies G(z) (inputs (z=Eᵢx) and outputs are both encoded) - internally called by `Sample` method
+    - `nothing`     : applies Dₒ∘G∘Eᵢ(x) (nothing is encoded) - most common for user interaction
+    - `"in"`        : applies Dₒ∘G(z) (the inputs are provided as encoded (z=Eᵢx))
+    - `"out"`       : applies G∘Eᵢ(x) (the outputs are returned as encoded)
+    - `"in_and_out"`: applies G(z) (inputs (z=Eᵢx) and outputs are both encoded) - internally called by `Sample` method
 - `add_obs_noise_cov`[=`false`]: When returning the prediction covariance, whether to add the observational noise
     - `false`: Only return the uncertainty given by the machine learning tool - most common for user emulator validation
     - `true` : Return the sum of emulator and observational uncertainty - internally called by `Sample` method
@@ -532,10 +532,10 @@ Makes a prediction using the ForwardMapWrapper on new inputs (each new inputs gi
 
 Keyword args
 - `encode` [=`nothing`]: For the output encoder `Eₒ`, and input decoder `Dᵢ` stored in the `ForwardMapWrapper`, we have provided the forward map `G` in the decoded space. Interpret the keyword as follows:
-    - nothing     : applies G(x) (nothing is encoded) - most common for user interaction
-    - "in"        : applies G∘Dᵢ(z) (the inputs are provided as encoded (x=Dᵢz))
-    - "out"       : applies Eₒ∘G(x) (the outputs are returned as encoded)
-    - "in_and_out": applies Eₒ∘G∘Dᵢ(z) (inputs (x=Dᵢz) and outputs are both encoded) - internally called by `Sample` method
+    - `nothing`     : applies G(x) (nothing is encoded) - most common for user interaction
+    - `"in"`        : applies G∘Dᵢ(z) (the inputs are provided as encoded (x=Dᵢz))
+    - `"out"`       : applies Eₒ∘G(x) (the outputs are returned as encoded)
+    - `"in_and_out"`: applies Eₒ∘G∘Dᵢ(z) (inputs (x=Dᵢz) and outputs are both encoded) - internally called by `Sample` method
 - `add_obs_noise_cov`[=`false`]: When returning the prediction covariance, whether to add the observational noise
     - `false`: Only return the uncertainty given by the machine learning tool - most common for user emulator validation
     - `true` : Return the sum of emulator and observational uncertainty - internally called by `Sample` method