From d73738c1c789751e31527f0097d5b16e46788ee2 Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Tue, 23 Mar 2021 08:42:58 -0700
Subject: [PATCH 01/10] first commit

---
 src/GaussianProcessEmulator.jl | 64 +++++++++++++++++++++++++++++-----
 1 file changed, 55 insertions(+), 9 deletions(-)

diff --git a/src/GaussianProcessEmulator.jl b/src/GaussianProcessEmulator.jl
index 00b656e56..2f855344d 100644
--- a/src/GaussianProcessEmulator.jl
+++ b/src/GaussianProcessEmulator.jl
@@ -103,7 +103,11 @@ function GaussianProcess(
     models = Any[]
     # Number of models (We are fitting one model per output dimension)
     N_models = output_dim
-    
+   
+
+    # TO DO: Standardize the data here
+    # Can use the full time median or some user define function?
+
     # Normalize the inputs if normalized==true 
     input_mean = reshape(mean(get_inputs(input_output_pairs), dims=2), :, 1) #column vector
     sqrt_inv_input_cov = nothing
@@ -115,6 +119,7 @@ function GaussianProcess(
         GPinputs = get_inputs(input_output_pairs)
     end
     
+    
     # Transform data if obs_noise_cov available (if obs_noise_cov==nothing, transformed_data is equal to data)
     transformed_data, decomposition = svd_transform(get_outputs(input_output_pairs), obs_noise_cov)
 
@@ -457,12 +462,32 @@ Note: If F::SVD is the factorization object, U, S, V and Vt can be obtained via
 F.U, F.S, F.V and F.Vt, such that A = U * Diagonal(S) * Vt. The singular values 
 in S are sorted in descending order.
 """
-
+# TO DO: Add truncate_svd as an input flag here
 function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, Nothing}) where {FT}
     if obs_noise_cov != nothing
-        decomposition = svd(obs_noise_cov)
-        sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
-        transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
+        # MFH, 3/22/21: Truncate the SVD as a form of regularization
+	if truncate_svd # this variable needs to be provided to this function
+            # Perform SVD
+            decomposition = svd(obs_noise_cov)
+            # Find cutoff
+            σ = decomposition.S
+            σ_cumsum = cumsum(σ) / sum(σ);
+            P_cutoff = 0.95;
+            ind = findall(x->x>P_cutoff, σ_cumsum); k = ind[1]
+            println("SVD truncated at k:")
+            println(k)
+            # Apply truncated SVD
+            n = size(obs_noise_cov)[1]
+            decomposition.S[k+1:n] = zeros(n-k)
+            sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(SVD_local.S))
+            transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
+            transformed_data = transformed_data[1:k, :];
+            transformed_data = convert(Matrix{Float64},transformed_data');
+	else
+            decomposition = svd(obs_noise_cov)
+            sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
+            transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
+        end
     else
         decomposition = nothing
         transformed_data = data
@@ -472,10 +497,30 @@ function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, No
 end
 
 function svd_transform(data::Vector{FT}, obs_noise_cov::Union{Array{FT, 2}, Nothing}) where {FT}
-    if obs_noise_cov != nothing
-        decomposition = svd(obs_noise_cov)
-        sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
-        transformed_data =  sqrt_singular_values_inv * decomposition.Vt * data
+     if obs_noise_cov != nothing
+        # MFH, 3/22/21: Truncate the SVD as a form of regularization
+	if truncate_svd # this variable needs to be provided to this function
+            # Perform SVD
+            decomposition = svd(obs_noise_cov)
+            # Find cutoff
+            σ = decomposition.S
+            σ_cumsum = cumsum(σ) / sum(σ);
+            P_cutoff = 0.95;
+            ind = findall(x->x>P_cutoff, σ_cumsum); k = ind[1]
+            println("SVD truncated at k:")
+            println(k)
+            # Apply truncated SVD
+            n = size(obs_noise_cov)[1]
+            decomposition.S[k+1:n] = zeros(n-k)
+            sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(SVD_local.S))
+            transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
+            transformed_data = transformed_data[1:k, :];
+            transformed_data = convert(Matrix{Float64},transformed_data');
+	else
+            decomposition = svd(obs_noise_cov)
+            sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
+            transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
+        end
     else
         decomposition = nothing
         transformed_data = data
@@ -483,6 +528,7 @@ function svd_transform(data::Vector{FT}, obs_noise_cov::Union{Array{FT, 2}, Noth
 
     return transformed_data, decomposition
 end
+
 """
 svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::{Array{FT, 2}, decomposition::SVD) where {FT}
 

From ba50b4c173776d904c14f5c27640def609bd9607 Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Tue, 23 Mar 2021 17:06:35 -0700
Subject: [PATCH 02/10] added trunc and standardization

---
 examples/Lorenz/Lorenz_example.jl        | 18 +++++--
 examples/Lorenz/Manifest.toml            |  6 +--
 examples/Lorenz/Project.toml             |  1 -
 src/GaussianProcessEmulator.jl           | 61 +++++++++++++++++++-----
 src/MarkovChainMonteCarlo.jl             | 27 +++++++++--
 test/GaussianProcessEmulator/runtests.jl | 20 +++++---
 test/MarkovChainMonteCarlo/runtests.jl   |  7 ++-
 7 files changed, 106 insertions(+), 34 deletions(-)

diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl
index c3f728617..7328e5935 100644
--- a/examples/Lorenz/Lorenz_example.jl
+++ b/examples/Lorenz/Lorenz_example.jl
@@ -246,9 +246,16 @@ end
 ###  Emulate: Gaussian Process Regression
 ###
 
+# Emulate-sample settings
+standardize = false
+truncate_svd = 1.0
+
 gppackage = GaussianProcessEmulator.GPJL()
 pred_type = GaussianProcessEmulator.YType()
 
+# Standardize the output data
+norm_factor = get_standardizing_factors(yt) 
+
 # Get training points from the EKP iteration number in the second input term  
 input_output_pairs = Utilities.get_training_points(ekiobj, N_iter)
 normalized = true
@@ -256,7 +263,10 @@ normalized = true
 # Default kernel
 gpobj = GaussianProcessEmulator.GaussianProcess(input_output_pairs, gppackage; 
 			 obs_noise_cov=Γy, normalized=normalized,
-                         noise_learn=false, prediction_type=pred_type)
+                         noise_learn=false, 
+			 truncate_svd=truncate_svd, standardize=standardize,
+			 prediction_type=pred_type,
+			 norm_factor=norm_factor)
 
 # Check how well the Gaussian Process regression predicts on the
 # true parameters
@@ -294,7 +304,8 @@ step = 0.1 # first guess
 max_iter = 2000
 yt_sample = truth_sample
 mcmc_test = MarkovChainMonteCarlo.MCMC(yt_sample, Γy, priors, step, u0, max_iter, 
-                         mcmc_alg, burnin, svdflag=true)
+                         mcmc_alg, burnin, norm_factor;
+                         svdflag=svd_flag, standardize=standardize, truncate_svd=truncate_svd)
 new_step = MarkovChainMonteCarlo.find_mcmc_step!(mcmc_test, gpobj, max_iter=max_iter)
 
 # Now begin the actual MCMC
@@ -305,7 +316,8 @@ burnin = 2000
 max_iter = 100000
 
 mcmc = MarkovChainMonteCarlo.MCMC(yt_sample, Γy, priors, new_step, u0, max_iter, 
-                    mcmc_alg, burnin, svdflag=true)
+                    mcmc_alg, burnin, norm_factor;
+                    svdflag=svd_flag, standardize=standardize, truncate_svd=truncate_svd)
 MarkovChainMonteCarlo.sample_posterior!(mcmc, gpobj, max_iter)
 
 println("Posterior size")
diff --git a/examples/Lorenz/Manifest.toml b/examples/Lorenz/Manifest.toml
index 97946c27c..8cabce0bb 100644
--- a/examples/Lorenz/Manifest.toml
+++ b/examples/Lorenz/Manifest.toml
@@ -1008,10 +1008,10 @@ uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
 version = "0.6.1"
 
 [[Rmath_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "d76185aa1f421306dec73c057aa384bad74188f0"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "1b7bf41258f6c5c9c31df8c1ba34c1fc88674957"
 uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
-version = "0.2.2+1"
+version = "0.2.2+2"
 
 [[Roots]]
 deps = ["Printf"]
diff --git a/examples/Lorenz/Project.toml b/examples/Lorenz/Project.toml
index e4fad32d5..f5c28e9b6 100644
--- a/examples/Lorenz/Project.toml
+++ b/examples/Lorenz/Project.toml
@@ -3,7 +3,6 @@ Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/src/GaussianProcessEmulator.jl b/src/GaussianProcessEmulator.jl
index 2f855344d..f8b6b4a85 100644
--- a/src/GaussianProcessEmulator.jl
+++ b/src/GaussianProcessEmulator.jl
@@ -21,7 +21,7 @@ export predict
 
 export GPJL, SKLJL
 export YType, FType
-export svd_transform, svd_reverse_transform_mean_cov
+export svd_transform, svd_reverse_transform_mean_cov, get_standardizing_factors
 
 """
     GaussianProcessesPackage
@@ -71,6 +71,8 @@ struct GaussianProcess{FT<:AbstractFloat, GPM}
     normalized::Bool
     "prediction type (`y` to predict the data, `f` to predict the latent function)"
     prediction_type::Union{Nothing, PredictionType}
+    "Standardization factors (characteristic values of the problem)"
+    norm_factors::Union{Array{FT}, Nothing}
 end
 
 
@@ -89,7 +91,10 @@ function GaussianProcess(
     obs_noise_cov::Union{Array{FT, 2}, Nothing}=nothing,
     normalized::Bool=true,
     noise_learn::Bool=true,
-    prediction_type::PredictionType=YType()) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject}
+    truncate_svd::FT=1.0,
+    standardize::Bool=false,
+    prediction_type::PredictionType=YType(),
+    norm_factor::Union{Array{FT,1}, Nothing}=nothing) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject}
 
     # Consistency checks
     input_dim, output_dim = size(input_output_pairs, 1)
@@ -107,6 +112,10 @@ function GaussianProcess(
 
     # TO DO: Standardize the data here
     # Can use the full time median or some user define function?
+    if standardize
+        #norm_factors = get_standarizing_factors() 
+	obs_noise_cov = obs_noise_cov ./ (norm_factor .* norm_factor')
+    end
 
     # Normalize the inputs if normalized==true 
     input_mean = reshape(mean(get_inputs(input_output_pairs), dims=2), :, 1) #column vector
@@ -121,7 +130,8 @@ function GaussianProcess(
     
     
     # Transform data if obs_noise_cov available (if obs_noise_cov==nothing, transformed_data is equal to data)
-    transformed_data, decomposition = svd_transform(get_outputs(input_output_pairs), obs_noise_cov)
+    transformed_data, decomposition = svd_transform(get_outputs(input_output_pairs), 
+						    obs_noise_cov, truncate_svd=truncate_svd)
 
     # Use a default kernel unless a kernel was supplied to GaussianProcess
     if GPkernel==nothing
@@ -193,7 +203,8 @@ function GaussianProcess(
                                                models,
                                                decomposition,
                                                normalized,
-                                               prediction_type)
+                                               prediction_type,
+					       norm_factor)
 end
 
 
@@ -211,7 +222,10 @@ function GaussianProcess(
     obs_noise_cov::Union{Array{FT, 2}, Nothing}=nothing,
     normalized::Bool=true,
     noise_learn::Bool=true,
-    prediction_type::PredictionType=YType()) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject}
+    truncate_svd::FT=1.0,
+    standardize::Bool=false,
+    prediction_type::PredictionType=YType(),
+    norm_factor::Union{Array{FT,1}, Nothing}=nothing) where {FT<:AbstractFloat, K<:Kernel, KPy<:PyObject}
 
     # Consistency checks
     input_dim, output_dim = size(input_output_pairs, 1)
@@ -226,6 +240,14 @@ function GaussianProcess(
     # Number of models (We are fitting one model per output dimension)
     N_models = output_dim
 
+    # TO DO: Standardize the data here
+    # Can use the full time median or some user define function?
+    if standardize
+        #norm_factors = get_standarizing_factors() 
+	obs_noise_cov = obs_noise_cov ./ (norm_factor .* norm_factor')
+    end
+
+    
     # Normalize the inputs if normalized==true
     # Note that contrary to the GaussianProcesses.jl (GPJL) GPE, the 
     # ScikitLearn (SKLJL) GaussianProcessRegressor requires inputs to be of 
@@ -241,7 +263,8 @@ function GaussianProcess(
 
     # Transform data if obs_noise_cov available (if obs_noise_cov==nothing, 
     # transformed_data is equal to data)
-    transformed_data, decomposition = svd_transform(get_outputs(input_output_pairs), obs_noise_cov)
+    transformed_data, decomposition = svd_transform(get_outputs(input_output_pairs), 
+						    obs_noise_cov, truncate_svd=truncate_svd)
 
     if GPkernel==nothing
         println("Using default squared exponential kernel, learning length scale and variance parameters")
@@ -302,7 +325,8 @@ function GaussianProcess(
                                                 models,
                                                 decomposition,
                                                 normalized,
-                                                YType())
+                                                YType(),
+						norm_factor)
 end
 
 
@@ -463,16 +487,17 @@ F.U, F.S, F.V and F.Vt, such that A = U * Diagonal(S) * Vt. The singular values
 in S are sorted in descending order.
 """
 # TO DO: Add truncate_svd as an input flag here
-function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, Nothing}) where {FT}
+function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, Nothing}; 
+		       truncate_svd::FT=1.0) where {FT}
     if obs_noise_cov != nothing
         # MFH, 3/22/21: Truncate the SVD as a form of regularization
-	if truncate_svd # this variable needs to be provided to this function
+	if truncate_svd<1.0 # this variable needs to be provided to this function
             # Perform SVD
             decomposition = svd(obs_noise_cov)
             # Find cutoff
             σ = decomposition.S
             σ_cumsum = cumsum(σ) / sum(σ);
-            P_cutoff = 0.95;
+            P_cutoff = truncate_svd;
             ind = findall(x->x>P_cutoff, σ_cumsum); k = ind[1]
             println("SVD truncated at k:")
             println(k)
@@ -496,16 +521,18 @@ function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, No
     return transformed_data, decomposition
 end
 
-function svd_transform(data::Vector{FT}, obs_noise_cov::Union{Array{FT, 2}, Nothing}) where {FT}
+function svd_transform(data::Vector{FT}, 
+		       obs_noise_cov::Union{Array{FT, 2}, Nothing};
+		       truncate_svd::FT=1.0) where {FT}
      if obs_noise_cov != nothing
         # MFH, 3/22/21: Truncate the SVD as a form of regularization
-	if truncate_svd # this variable needs to be provided to this function
+	if truncate_svd<1.0 # this variable needs to be provided to this function
             # Perform SVD
             decomposition = svd(obs_noise_cov)
             # Find cutoff
             σ = decomposition.S
             σ_cumsum = cumsum(σ) / sum(σ);
-            P_cutoff = 0.95;
+            P_cutoff = truncate_svd;
             ind = findall(x->x>P_cutoff, σ_cumsum); k = ind[1]
             println("SVD truncated at k:")
             println(k)
@@ -562,4 +589,12 @@ function svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::Array{FT, 2}, dec
     return transformed_μ, transformed_σ2
 end
 
+function get_standardizing_factors(data::Array{FT,2}) where {FT}
+    # Input: data size: N_data x N_ensembles
+    # Ensemble median of the data
+    norm_factor = median(data,dims=2)
+    return norm_factor
+end
+
+
 end 
diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl
index 14d8c5b0e..56b2794d0 100644
--- a/src/MarkovChainMonteCarlo.jl
+++ b/src/MarkovChainMonteCarlo.jl
@@ -74,16 +74,26 @@ function MCMC(
     param_init::Vector{FT},
     max_iter::IT,
     algtype::String,
-    burnin::IT;
-    svdflag=true) where {FT<:AbstractFloat, IT<:Int}
+    burnin::IT,
+    norm_factors::Union{Array{FT, 1}, Nothing};
+    svdflag=true,
+    standardize=false,
+    truncate_svd=1.0) where {FT<:AbstractFloat, IT<:Int}
 
     
     param_init_copy = deepcopy(param_init)
     
+    # Standardize MCMC input?
+    if standardize
+        obs_sample = obs_sample ./ norm_factors
+	cov_norm_factor = norm_factor .* norm_factor;
+	obs_noise_cov = obs_noise_cov ./ cov_norm_factor;
+    end
+
     # We need to transform obs_sample into the correct space 
     if svdflag
         println("Applying SVD to decorrelating outputs, if not required set svdflag=false")
-        obs_sample, unused = svd_transform(obs_sample, obs_noise_cov)
+        obs_sample, unused = svd_transform(obs_sample, obs_noise_cov; truncate_svd=truncate_svd)
     else
         println("Assuming independent outputs.")
     end
@@ -173,8 +183,15 @@ function log_likelihood(mcmc::MCMC{FT},
         diff = g - mcmc.obs_sample
         log_rho[1] = -FT(0.5) * diff' * (mcmc.obs_noise_cov \ diff)
     else
-        log_gpfidelity = -FT(0.5) * log(det(Diagonal(gvar))) # = -0.5 * sum(log.(gvar))
-        diff = g - mcmc.obs_sample
+	# det(log(Γ))
+	# Ill-posed numerically for ill-conditioned covariance matrices with det≈0
+        #log_gpfidelity = -FT(0.5) * log(det(Diagonal(gvar))) # = -0.5 * sum(log.(gvar))
+	# Well-posed numerically for ill-conditioned covariance matrices with det≈0
+	full_cov = Diagonal(gvar)
+	eigs = eigvals(full_cov)
+	log_gpfidelity = -FT(0.5) * sum(log.(eigs))
+	# Combine got log_rho
+	diff = g - mcmc.obs_sample
         log_rho[1] = -FT(0.5) * diff' * (Diagonal(gvar) \ diff) + log_gpfidelity
     end
     return log_rho[1]
diff --git a/test/GaussianProcessEmulator/runtests.jl b/test/GaussianProcessEmulator/runtests.jl
index 11622af28..64ae66bfc 100644
--- a/test/GaussianProcessEmulator/runtests.jl
+++ b/test/GaussianProcessEmulator/runtests.jl
@@ -49,7 +49,8 @@ using CalibrateEmulateSample.DataStorage
    
     gp1 = GaussianProcess(iopairs, gppackage; GPkernel=GPkernel, obs_noise_cov=nothing, 
                    normalized=false, noise_learn=true, 
-                   prediction_type=pred_type)
+		   truncate_svd=1.0, standardize=false,
+                   prediction_type=pred_type, norm_factor=nothing)
     μ1, σ1² = GaussianProcessEmulator.predict(gp1, new_inputs)
     
     @test gp1.input_output_pairs == iopairs
@@ -63,7 +64,9 @@ using CalibrateEmulateSample.DataStorage
     # Check if normalization works
     gp1_norm = GaussianProcess(iopairs, gppackage; GPkernel=GPkernel, 
                         obs_noise_cov=nothing, normalized=true, 
-                        noise_learn=true, prediction_type=pred_type)
+                        noise_learn=true, truncate_svd=1.0,
+			standardize=false, prediction_type=pred_type,
+			norm_factor=nothing)
     @test gp1_norm.sqrt_inv_input_cov ≈ [sqrt(1.0 / Statistics.var(x))] atol=1e-4
 
 
@@ -71,7 +74,8 @@ using CalibrateEmulateSample.DataStorage
     pred_type = FType()
     gp2 = GaussianProcess(iopairs, gppackage; GPkernel=GPkernel, obs_noise_cov=nothing, 
                    normalized=false, noise_learn=true, 
-                   prediction_type=pred_type)
+		   truncate_svd=1.0, standardize=false,
+                   prediction_type=pred_type, norm_factor=nothing)
     μ2, σ2² = GaussianProcessEmulator.predict(gp2, new_inputs)
     # predict_y and predict_f should give the same mean
     @test μ2 ≈ μ1 atol=1e-6
@@ -87,7 +91,8 @@ using CalibrateEmulateSample.DataStorage
 
     gp3 = GaussianProcess(iopairs, gppackage; GPkernel=GPkernel, obs_noise_cov=nothing, 
                    normalized=false, noise_learn=true, 
-                   prediction_type=pred_type)
+		   truncate_svd=1.0, standardize=false,
+                   prediction_type=pred_type, norm_factor=nothing)
    
     μ3, σ3² = GaussianProcessEmulator.predict(gp3, new_inputs)
     @test vec(μ3) ≈ [0.0, 1.0, 0.0, -1.0, 0.0] atol=0.3
@@ -127,15 +132,16 @@ using CalibrateEmulateSample.DataStorage
     @test get_inputs(iopairs2) == X
     @test get_outputs(iopairs2) == Y
 
-    transformed_Y, decomposition = svd_transform(Y, Σ)
+    transformed_Y, decomposition = svd_transform(Y, Σ, truncate_svd=1.0)
     @test size(transformed_Y) == size(Y)
-    transformed_Y, decomposition = svd_transform(Y[:, 1], Σ)
+    transformed_Y, decomposition = svd_transform(Y[:, 1], Σ, truncate_svd=1.0)
     @test size(transformed_Y) == size(Y[:, 1])
 
     
     gp4 = GaussianProcess(iopairs2, gppackage, GPkernel=nothing, obs_noise_cov=Σ, 
                           normalized=true, noise_learn=true, 
-                          prediction_type=pred_type)
+		          truncate_svd=1.0, standardize=false,
+                          prediction_type=pred_type, norm_factor=nothing)
 
     new_inputs = zeros(2, 4)
     new_inputs[:, 2] = [π/2, π]
diff --git a/test/MarkovChainMonteCarlo/runtests.jl b/test/MarkovChainMonteCarlo/runtests.jl
index eb6efd280..a9d66f245 100644
--- a/test/MarkovChainMonteCarlo/runtests.jl
+++ b/test/MarkovChainMonteCarlo/runtests.jl
@@ -58,8 +58,10 @@ using CalibrateEmulateSample.DataStorage
     burnin = 0
     step = 0.5 # first guess
     max_iter = 5000
+    norm_factors=nothing
     mcmc_test = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, 
-                        mcmc_alg, burnin)
+                        mcmc_alg, burnin, norm_factors; svdflag=true, standardize=false,
+			truncate_svd=1.0)
     new_step = find_mcmc_step!(mcmc_test, gp)
 
     # reset parameters 
@@ -68,7 +70,8 @@ using CalibrateEmulateSample.DataStorage
 
     # Now begin the actual MCMC
     mcmc = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, 
-                   mcmc_alg, burnin, svdflag=true)
+                   mcmc_alg, burnin, norm_factors; svdflag=true, standardize=false,
+                   truncate_svd=1.0)
     sample_posterior!(mcmc, gp, max_iter)
     posterior_distribution = get_posterior(mcmc)      
     #post_mean = mean(posterior, dims=1)[1]

From d652ab69b57ae65d9ceb6cb8ba367fad9245cced Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Tue, 23 Mar 2021 18:26:28 -0700
Subject: [PATCH 03/10] working lorenz, need to debug the truncation

---
 examples/Lorenz/Lorenz_example.jl             |  13 ++++++++++---
 examples/Lorenz/Project.toml                  |   1 +
 examples/Lorenz/outputinput_output_pairs.jld2 | Bin 0 -> 17704 bytes
 3 files changed, 11 insertions(+), 3 deletions(-)
 create mode 100644 examples/Lorenz/outputinput_output_pairs.jld2

diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl
index 7328e5935..524bc1601 100644
--- a/examples/Lorenz/Lorenz_example.jl
+++ b/examples/Lorenz/Lorenz_example.jl
@@ -89,7 +89,8 @@ end
 #logmean_A, logstd_A = logmean_and_logstd(A_true, 0.2*A_true)
 
 if dynamics == 2
-	prior_means = [F_true+1.0, A_true+0.5]
+	#prior_means = [F_true+0.5, A_true+0.5]
+	prior_means = [F_true, A_true]
 	prior_stds = [2.0, 0.5*A_true]
 	d1 = Parameterized(Normal(prior_means[1], prior_stds[1]))
 	d2 = Parameterized(Normal(prior_means[2], prior_stds[2]))
@@ -219,6 +220,7 @@ ekiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_params,
 # EKI iterations
 println("EKP inversion error:")
 err = zeros(N_iter) 
+err_params = zeros(N_iter) 
 for i in 1:N_iter
     if log_normal==false
         params_i = get_u_final(ekiobj)
@@ -227,8 +229,10 @@ for i in 1:N_iter
     end
     g_ens = GModel.run_G_ensemble(params_i, lorenz_settings_G)
     EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens) 
-    err[i] = get_error(ekiobj)[end] #mean((params_true - mean(params_i,dims=2)).^2)
-    println("Iteration: "*string(i)*", Error: "*string(err[i]))
+    err[i] = get_error(ekiobj)[end] 
+    err_params[i] = mean((params_true - mean(params_i,dims=2)).^2)
+    println("Iteration: "*string(i)*", Error (data): "*string(err[i]))
+    println("Iteration: "*string(i)*", Error (params): "*string(err_params[i]))
 end
 
 # EKI results: Has the ensemble collapsed toward the truth?
@@ -255,6 +259,8 @@ pred_type = GaussianProcessEmulator.YType()
 
 # Standardize the output data
 norm_factor = get_standardizing_factors(yt) 
+println(size(norm_factor))
+norm_factor = vcat(norm_factor...)
 
 # Get training points from the EKP iteration number in the second input term  
 input_output_pairs = Utilities.get_training_points(ekiobj, N_iter)
@@ -297,6 +303,7 @@ println("initial parameters: ", u0)
 
 # MCMC parameters    
 mcmc_alg = "rwm" # random walk Metropolis
+svd_flag = true
 
 # First let's run a short chain to determine a good step size
 burnin = 0
diff --git a/examples/Lorenz/Project.toml b/examples/Lorenz/Project.toml
index 54a9e6253..c3b779a9b 100644
--- a/examples/Lorenz/Project.toml
+++ b/examples/Lorenz/Project.toml
@@ -3,6 +3,7 @@ Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/examples/Lorenz/outputinput_output_pairs.jld2 b/examples/Lorenz/outputinput_output_pairs.jld2
new file mode 100644
index 0000000000000000000000000000000000000000..546e4eb36205cd9cabf04adee690da8059fc90e6
GIT binary patch
literal 17704
zcmeIac{r5s|2I6$*thIiX2~{FW-K#<ki4vERVpRQk}cT^rL<^KDwV`1Nk!3WDN&TP
z$r7ows7TTxl%!B7&*{_Wd*9FVJMQ~;+|QrSb9~O@7}s^pbzbN4KHtmh{W`DV;^=I~
z^Ye`g@Mq1430>jKa`cV%WjRNz4)cxXvepEwj#?2B&N4AFH8Q2JIDe{`8krdhS>{5X
z-->9~ET@T-f6M=8C9uoI(OFJfQ(9C~l0p&vGmX)I3diZsb&36-qw;#&CRZ0nH&Mz?
z>QoBMNp%*v6Yd)pK%rPFicrRpD--e&As-4wk`n2=n%s*fN<<lFNxm7oZ}ImX3PskA
z^d}-ik*3fnA!C~79lb6xfI<}=YhbL54MvW9YQ`=|`P1m1^Z)!c&ybq+=X{!j>#V<P
zjP*f^$UpVRS{N%15vv1?#&mnZSTAF(|7p+k@4FO=Q+QNBm|tkXOyAHj-|)GsBbEn5
zMFm9t{mA_2h}FJ9e;#*;2#@w%5gxGmpXaFmv?Ke^_22zCher5Dn+wOjqWyE9NcZnD
zUGdLR?4R$)#4%R>Px_zhyJB^K|JWq{=P1PgH|_svl|t#2|IgXK$A2a8uLS;;z`qjs
zR|5Y^;9m*+D}jF{@UH~^mB7Ce_*Vk||C7KUo+6M+Cp-no^vBZgD3Znp_<Q@4{Q3VS
z4KYR%{D+1Zqb0_uga2I9$ocsEsjL2U=qN_+Q>8@6w__El(iH0dWHco8#`pjhWsEEM
z&*M~b{$Goob&2p#b1l^Qe|VNLR%1nYWK47vh4SCdCn6^L?^=@oZLR+!@@0&g_|GuK
z$oc<$*5c$^cFNc!#}V1eQUX^5g#Hoj7>g72KUf}7!o!RSVyry>hno4%QQX0F+n;F)
z|5s8=|A#4&R2h@nI7)@mpJvYJ{5h(dr@Z<5agkPqf1W2bMG3Jp^ryn#D=>CSRiTKO
zvTXlW0ty8bp$tVhN(3*Nqd;2DQaJ2+`G=v84uVc71ob_a#F`t1mf<RDkPPT}&s3Fx
zO0GuuBDy{lowit4^+_OW+`FG_cOA?UqpC|84%>$AoSoXnWaFi1Zs|jcG)`qH7SNjH
zAaLd6Pvpy@({`1@`?qvt8PBrpY8<ltk#LvyNo@jn5sq5512Sl@%yV+BSH<u_JqwfK
zueM@YPErk5b<iHjrAqT8P<gnddAFkuKHT2->9CqIu3wo|{@Z>6%G<Vy#%-6u<<AQy
z2A*Pp@@-OSPLmScZtP4Luw+4YU{mX)^>PSl5AvRLo(|saPi6{C8Z@1R(JzN|!FpO}
z)bl_Js}i*@YHiR4o*xYVm_2H{XQsd9^#V=2yBz2fK~=;DwWT#%9?;Ra(WN2niYy{K
zGka>MvcV`)ZrfR`h^su`qRfx_P!;+3)9$tydZx~K+Ix|Wl1h~ux6MUivngefpMy5|
zF&j!;Xc8z-`#pZ$8Wx<gzX!j*C4-zTzgnJ&Xu^a3%81e=fxGRB^&vtP1UIdJU#=$s
zeX+X7!89Eh*zGN$71NMW@y-3vD;D1PJIT#DDS`f)vrS)*Xdu3K*5|p&vPhaz;MFQK
z9!_eP69cu!Ay`HFbgVB6?nW`^3{T49n#K3M>M?Yj@*f)PpE_!LPc_Aa8l-`Y)p3(+
z?Bx-FV*1K71FSJjR}D!TuvK~JDSdsV4kn5`@F@F2L+o4i2+(!$I5ufUo{St8f6OY)
z`@zJgC^4~*6D6RzHAO_ckO_C0g~xZ@SBBBfV}+d?RZ;NH=jY5uN!)rAa8+%mHtez%
zeqp>5MfkeNEh{?|Q6s14mQ_Z<$*cC1*{LcpUugSsU9UXutrAt(w^J9QAvRNZmu29c
z6mzRIhyiiW?71%rM{RdcS<0DvQU|{WyAH;>P@#Elvx}sn4!S>X8#vo6g&jq$H#`c}
zkTY&+Nbf8~1P;fp*%YCVr0Wm+vs*<lu=lF!hc*_N`XLu>w}_#^x%Ov^yEa7byS&=|
zRRUAgGJRjBG0_rmAW%u4hV++_{mlU^oS1I6<7lD;sxDm=-*$;XuI2if$>Io<$$vUU
zk%eU^51%jh7lWJX&|g8vm>50XwRum!IINOoAIzSp4W~kd<bfV3ths%~M`<Y&1?xt~
z7Y~b~<Kx!#ho-QgpSLT-mqA4it;c@iVkSPEEYjUGQ3Ne3-^B2FEclGL4lH~qhJ_z(
zY;>M6amygJPL5o!S$O%ZQa2Oh@*R(sU8KU<Pn<SKhO`%#_HvGn7?uaR6>h!DgyvGm
zktL6ayx3E8HncF&JQVcu;BryCoWU=-smnyqPlKdO0vg8sT4NPo%S1v&8C67H6z$<9
zZ?fZ=s0fIQ`*Bko>DH2m)j2HKT-ouL<#sW&rMIcS*vi1~RVg2?v8Xt)U4wVAh=FOd
zen@0xQ<28)8F;vsff0@{FE&67p*n5b1H@RE`#h+!+Cl`2YG*5TB{N|sd-9>@87f3x
zSv&}y&BVcDi<4M6R5%`XQmJM$ku7jAbMB{Me8)HYL0uMdj%>Vpc!Y)<gS0Il?lG}F
z`mC(;VNvX<m1+I!EE72;>34EEY1pSf+}>Tm#C$*JoT}4esNJG8&OesuD?6jqO^(dh
z{iRC$UMALQUEKZ6hKi%xEtP`IS(sCDNouvR7?#Z4xPMm?6P2lZ;)?G{qs~xTvw9i}
z36F}GxStfq^5FO>Kgsy>jjPsN^%uu{ksT~kQx={V-J5avHx)m9CX@`UC*wP#wyJfG
zIN0@R#hzL${3^fdkQOC@cOO~v8YwJDpITed_*E1a4jfO&Ol08Yl%e&#>SB2HB~Le4
zlZ77>YBa)#fBhtN?m?0~3vcfy?%i}&0vZRTm<|t_I6ic9@_9EAeD>39q5G2c5-+G;
z>MDlhbeC+9^|Csj>Cvhrjx&W1?_`nwqwkHc3^0?%G8SWM{sI<;j_Svhc8jB1ZnJVT
zsV|$c(V~3?6~2o#CSS5(LH=4>aYr{5t?VhYEngD*ANiR0)slw03McJC4O!4unj~$q
zjfOK%RM+(sF)<^-G?*ixVtm}Ph+Sm8DlJCG<uy|=dDcY_K@}5PuMYMVvB~<L9r;%A
zhUg!6;}Kh;|I@w?O^ct;#ExFuik9=D7^pb4VzGdQSpf!#lJ}?>U_Ea3A7r8`Qa0$D
zkO~R+Vk-uVh13rPEGcUmuH2lUH>|;e>5;GpZ%)t<wdHDxn+OYkwe3<Gy-!8cledN<
zZA>UJkG%gSMT5*9k?dtdOsHQ8Pl;mDkd-(5em0GTNf#bp6r<1}csMzFehm|wW`?Yd
znNG#q-_jc`$$AO2nmegEH29R?G^`h8A@-NgJv|2+9#waj_7izLF>82hc$fxTZ!RO8
z=*JK1z_=GzMX~+cbQutTVQW1zP`ptTso}A+9V(ef+i;EdPKgHIG>xfKiF{qpwHa*A
z5W&OL-=a3ep2O|$E)F90)N+U0K9tSGedh5(w_BpnwSS~2p38(8mHnpu1QoaQ9FFYJ
zW5Kvx`^K?&8g^Yekvf^!yAeO{tA_y<9ZNYUS1Gfga(iLY<K;AbovPm`X2QbYyV(1B
zIy3~{Rh42i5c?64o*Oik2FdzwFU!iA=%3|SM?XfzflGl)XDno4mjElSYS6$`3v(Y-
zVnMC+M}VuM7!u_4j`))PX+yMQ&YCnNp0?>;Nv>bp{Z3b?BmT$pX^rOyk!Nvpy>Sn@
zU%Vn&sgBHt_AWbn6&1PTb-WzdEIiy1H#xqF21fhW0E!$7j2@#42Ir}GJ)f1+bb*Px
zUn;(t92Z5oOx0BuvB%Q5TAA_ZX_#Bd@M^PW!Oi9ELT?wM|But}{37;$%|=RZD4z<;
z{fsX$rA++JU!O5`7Zn*|6;meHGx5`DepOlo(Wi=iqd5~;xZLwgr+OL<PaVZ9y2yGK
z*VLNpFQGyGOp}ui(Fb?WkPMF`qE9Pth#sJ`FdTa3QK%{v&(Bjk%7}kfiPDnJ?5D!+
zxrKg_Gz%B5sQB-p(J=kSg0@FQ-~5wpCajYp`h2KKW(Ajp58`S1^bsm%WCY1iAo85z
zPC4jvl?p5WhcYA0OzfF-@#{W0DhiJ>Y-M*b;Um80@VWz{h}P`vd#l4jRdrfyjUo+4
z_ii59OZu~Te&WPJZ5qtoJEZ1lv*2g)PT<0&VWZM=n;ueMGf3SbA(00Eb!!te$@tYY
z4TnY!Qo&c<P*6hjxh03^cf5}ZSMTbdQN%v4WDWfA45T7Y;`NTJ5-hxWvN(6!Q!2K<
zIzO}b8WS-StW@Rw$b3V3ulT)ZLUw(ptidW84)dE=f%q$K<>_m+o2a-jdU4fwD;D0R
zaAF13G{{vgqs}LIVmhTS?fsxAG~`oq-;N{k#qMrfx(*dX=c@7K2AMx~GN<Yz729^s
z8_oX6#H89e8IA%`eATM|($vqy*$?v)|5`)juOQ;{>NXSU&oAd}BL0W|s=+C?k%`6$
zObHh;8mx_CGWo>+?K;@tUQH+Q!#nwb#8)Qv4K*hgD^hW4^GTUPDHdjA$?6&JqhZ1)
zr9&Y^AEFN!e(xEiBH1|o@!R`M>^gm_CzHevg(Uf@ETTVGUd`HO+CT;E<%@y0#w^^B
znz{MiS1QDsa?HPx`52!}Q5-r>!yNZZ4}6II@7SAK+*(P6(`-?WB=Jw>ABCrk-H5;S
zmS`&?{>ArV*rEp^1iv<I`7oV~&w1k4V`UvwQ0!mN%$>-B>aJh2g(WoXe{*v*iui}g
z8`jTB;}g6%yYXfb!5`BL#J$ozMUkVwp4~i$h3<K-pDs+M;dSF=`CpS+@CsSdcV?9s
z^v|Tx>@PAgUM(g65Q+b7XWq=&w~q;(`uDeVdPH${vB$F(vi_{sT}AhfP%(Ywb+Zg&
z4+m#Ae79LfLvMBI7Y!8_UWPu(7`;VBMRdg6tNBE~w<g}2J&uYsOZwi7>DTi!$+Q#1
zf87t+?3_#HSF}=HL4(M@H^SSbjKRXTis99dX3|h+6xo^IMB>8?+vVQGo{srk*%n3g
zb-u%@c|2npw#BUDn-Kl$?Oi{+_YM_QDyPVu#DlT{E93KVRIGP;(w62<@Ymf-=~LXv
z_!kXX+Y@=mDCe7`GHCGr>xZujkymC_Ps<{rpV|J+*3!g3(MxB4m{U$gVrk&UOoF#8
zZr9oi>}mKb+eh~Ru?MeIu@~bfQ&GP>zv&6FU-9i;qo0Yrz4KP`5V^xdWWtNZZ}q7t
zV|Ja1Bl@11MtiZ0^tbU*z3m&KA8kvOb&_JKSj3jsGkMFz?E42_dl7$fzB_kK7r{e|
zr|wFvCjRY}e46C~M}j|owsm%Xqv1$%w^lP*e_b(=Wy}#`9}Z>?u4MkpXWA^<K9kr}
zz__mTS4{YZbO`Sf|KTrNuRToU;ZYfFyzPM~5+~Bcl6y#hA=`H3XOZ@2{G5E{Gx47V
z&X*r_Q4t&F6DmvQ_geJ#ok9|yQ$zh@gB3`;o#wnR<~bFSx9fEV$a)>N`29R}R1~{w
z91ispdHRYO)P~j5FiGmor8=Syj;@}<U#c{0a;wz&LgL5l&eO_Ui9g+ulPD!i^sCn)
zM(Y>B6ExaZr?Zn-Fkmgxx3;0dvm;ZoM4W}KDI9xx1P$C=t;_CeER?L?^SY;p3Ljtj
z9s#*N9ap@WdY20AiY1IPQvXR!(SrTTG@Ogi(@rAund@P{{p&d@l6*Ck$J?+lUi7+F
z(P<ir4rd*0A^xbGKka}j!3Q0|n+Bhe_!+EuuQ=$kD2BW!(_8o~WTbvt*K9yT+~sk(
zgTy~wa>xi#jTD3FLZ=4bR3-+BzgI+?(eO;<PWWDOo%|{}cA7dBzc0k-`+gwyQLJ{S
zllZ^m%=uHh$owRw92bX?_HI@8x~`Wc{xoo*^w2gEZ({Tw3J5+`klWXEmrupJ3C~7m
zlKMXB{j3LzY3MTZqW{VCH!;OT|K?ZzOQt_Y!i<&wZKm(YB$+;#UgE}28unV@ch~h`
z3wJ&wO@2jdT+l&_l=z_mCJ&Xd!mxXnb+NGj>=}G#VnX)LmP(f)+s)qR>1}pOsMxS-
zse88#mK;jCcyp&YMo(PI8(cF9dN0<jUbaL)H2?QtjlLOXs9fx^C^toyLrYs~AqQ!P
zdv1$eF+xS@+9$cofrI;RJX-R^5;|JR=YsAS;c%7Q;8RmGd@*iX;^n4`3RREj*JqR=
zo4Yv(AFkSR^|lS1N~XZDf1|aFoC%tr8bv0Tagch4vf7hqV&Hpa4^Eyw{HLGHwJV;C
zC;8q}%eI(6zNPaqOHUIE9@c(YpvFPGy}nDrO&zcrq`pTFYT=i@T2t)-B^)p{mt8$?
zJeJO-{|-HFi7LOCIWw&Jc-Vcsxa+nqT1$h2?0b!{DNQ$;KV*aLwm#iepZGA8xjyjv
zg&O85ERZYwZipN&-by!BONiy%J2_fwf(iP$36(0^crv|4+&yt3;+>a9&`w&xI&K}W
z%Ebf{kteQsY#)a~dy$&zvj*5-wq$l$h6R>+8CZD@7-6sGxp2=6RcPG0G_2vn#(O*6
zp;hzw=nQMQA!j%VcT=-e{cns%NXDay_lI<$cOcDX^Xn-nop7`NkeM0YeTvN1Gc`cg
z<{$iOhAD(28+KW2GzM+z$4ZY@HmLHJCT*`Y5IQS4b>eGXgngQ~+ON+P9M!$%CBAIT
zu~d|fqlsYdm3dq1-ZG%bxZD`gAVim8+v>#M6EXkr`Ce(EA|`$|nU%Ck7rXD@`U{2*
zxGE{e^KiC?dv{2u<zWu2Q@>D7NSmQ7Th?r}dMc!f4Ej7LPsWA#r86V+c#v`BUQb^q
z#1Y4PnwmC%_L>ZrtXymC4j5MyZOg}UiLd-jM+*#`jZf%2WQv2aHS<l6a`5WXaJPrD
zGU~Q1_c)|W7%9$B)09aT$hi`BNBW~N<X_jc?ex>Yvy^RvIocfPu6m#DwA2JCGfRU<
zH|ZmR=hy4FO%_w;zxO)8A{q3$E8<n3^zq7?>!xL?jWm01d(VR(w!t&acbl7y!;H3k
z?;Exj7-`cyzMpG~y6yMNEsy8|hh)CUL~(Kc!l%ppr{<t1?u<>+Gr`=e>u!E)q+{)2
z=E0-4Igr<iWUm#N!0VNa+Og^Skg>5%$R-x9{nO#5Uw}4bU5|aZvc(kdwc-p@z4Wl`
zgWL*J{Ing-uJn%bQ$e9ySo^6Y3wX7&57kr{!%jroCUmws-kHXB_3Cq=T<)73aLEKW
z0w#p+PvhX&cI&Eac{)0+${fv|^fCO+p3}XA4_|@A%x^S(Y_1DlxMi&<n$;Jbdg84H
z{%+CshQF*a`*X^P{O{&?(P(@1);Iy^8PhjybTvnArD}TpPcr`-Ki(^8n!;{>U5J*p
zKC1hdoLV>A2#1E}Jum3D1}~*!Jd-0pf6kI?qsbGnx&6v-9~~1Ey^c$s8Y9G@J}0oQ
zj|=}9O+g#dHBqizQqWgv0RM(57cvH{AhXRe?9z231O?=_ZXJ=t<c9&>x)yp+f4)#{
zgQ6a$N9}k^-^oBnLES*&kr7)z%QI6QDi!g3@s_TG`%K`ru$4yH#D*-FopbG*I$;{5
z+P2ToMbE*oyYp%^usQzL18rwT^qpVQ=u|pn%TQ~7)>tcp3p1x&T<ocfxf;frcOGj(
zMCXu#NLhz%g4HE;&vpg$df;O@rodNMxysQ_fOYFDl#b}?!tCw+hqF6OuoH)TALy9E
zu}pA(rTHYt{-}|fTc?dnGi~H@nS@cxZS=YwLF6~n{M-ozKK52kew;Gb04u+l6!D^s
zG3(_<F&CT32xK)}@mXn#6se8Zn)6ubizwd^!s226(~!pc%acfktFY#gmkw@ENK(7l
z(``FNL)9Adl#poX#*cetjT6c7M_)yl;n>ulM=TS0IC7IN`BhVZ$wk(#Lu<?+@jGp;
z@i=3g>SN~D+mXyQXS#1#7tznr0{8rIKAcr5=bW9!#$8_U*69t3@HE`^mSd;`hGz*)
zC4|iT=epf@E*fI+qv`sRdODghdwlykHm3Bgj{f-73{zgjhfVy&gW9;!m!{U*IFVs`
zO2V82*{oZIdm;spVb!&niyIOq<oZtP1Pw%<-CrakYJh~(<8#{|2r*n?a3uUME=*db
zn^Uh5JAAt~T_lnXTK2pbWgE?K<4Id~l?fM5)Cz|0Zr4PytZ~U;o&nCT?GkS~z=yR~
zjNifcY=jyG*_TTwU}oH4=xd5LG^opM0t)!(8_=yVn#%=iqP%{{GdePjlJ?DeZ-Dy5
zxsK6fep{2@PTn|~hjy!zb@Tfe*z>N*=gvtEPDv%6@!ZA7$Cf1p#cvI<dSKvK<$MO_
z#s`+0br~S!>6?l$iWzd&)b}f8^04!{rcN-Aj%h}}oQK`mu-AAgdE8Nm76bW~7pJ+<
z+5bv958BwC6R)=A1{+7>l*D%&6Cmt?yv?})E;e|(&VMk%K>gYE3%u5{;rv5AY@eDC
zPv1Q|uql!U*-_OLK`R}yMx~kYhuBysQu^_-g#fcNJFm@h;Nn_;rp(C(I{5QiZ|jZY
zz~#NyBG!FASSmEOvNspsUGFYmHk*lw&T<9W)*SR2q#GR95un-dw^pgBAto0}-h8L3
z4Tp{O1uN4H(6ZPlp)!GwW`~><>sN-TyXH4Rlpv?V(p6~zyV=-94J;p^2~pJfS!73y
zA+(&oAGc4{#J0&FBX1NN;K@d+=PwH({CB>Z9~I0)tE0G%+^9D4CA4B{j5+Z8^6A``
zWIl|y?uhnqnuI)jI504*g)1vB7lpR!!*t2I=;fbHabxwz)8}#~Az;US`-731ct3ft
zCY#(pQFv$j_89_LhDkhg51#~)mlZ!gdTHSChSf@@jrs_4;U-2J3$Tc~H%(28hg8Wq
zN972!L`ybIY<tVb-n1<TO<MTKyEa~7e!n3errg<?&Zpy<p3?06tp-?k(Pq)&cmcFc
zeg{0+z=iw+ho43RS~zom?eX7wY_uo2?Y*NXM5xr6y1_0(Y-fAzlf12ki476O3X*Kp
zJa;iY|4;zcNkRKo5WJNttyF$gnt`Vxb%)FyIM5Wc3b=18!0i>EoTqzm5pNXU=9s~N
z<X-1Y-&w@I^0Vj9n=8PzvuhGqt%j&yCoSDqr3Iy^C6WEd*tk)ZwO9PE0G(g5<#*<A
zF~QC1tA)8XuKVx0@x|W&4tIicntTLUl9aV1R-B7Z<t+J>csfouh_71Z&c-O=M+3U~
zI8DhKzilNK+NqD5?9%D5l?pKQQsLl-jl0qcfe_nUTfRSR;=*b8?E06{+BoSvS9Y2y
z8}Ih-Ug0yt48y*Ci#Jtskxo+(y>d+rdhb`y&sJiiHhg?Tw2>KVFC0i|Si^(2+Y6bb
zx!Q1ysv24<!@>6z&dx$pA(Vzf0z%ph34`S#F@d9rqiKag8zR4Q-&d0#OADc`N6BiV
z^3YhZXvccO+*~>n9-H%lgZ`C_tR?<JXc+59bi{F?aYKAk{bx<wYk%Et@6AS}dFi3u
zvSxte)*oWmd00_X>Kdh^h4P!HFW%hG#uJB$@8XXN@U&P%p(fK1*LKwyb@wn3{Dt}9
zXek?O>t{aEc*V!y(3=g>v6C=??lXVxat#DLaKKH1C$tZ__;WhVVb<I6P9&Fyvd>{I
z!&Pd~U2eon4l;m_^r@P?d1QQ-8pk*K@}T;APo!d?CPIdMUcbJ_#>Vq<PZtF8v1ys>
z*UK?n^!uIr{Wgy<3;Iqr!Ui@T?>es6Li9T+`^ZauV=jJ_*{$UVXu)k=pw%U^zV)Xa
zca=>sLw~<jz3l=cNLfeLMh$Dj`BioId{qu=uADidK`{qq$4_fhbz?Lt+}vHBJsz7c
zOp>YeV<XdMhUrMP5bLgfk#MQtVdT`kzKT!U*l8ix%pv~Uf2Qm4kTU`_evsfRxAI_8
zRM)?1923zx7uP$U<=}*MV%oT2Av}%7e{)pf!G1;giZ#Y$r=~w*lhXzc=5Mx{tC1!^
zrBwsda6S)??Z%Z^%jnqKr=5{>g$?etJq_W8W-w&fkM_UkVr%V|M}j;CTGY>Myh_&h
zhpa)~xduLVE>~vlHsyl5>IqXuO&hzloDBAp=hts3UV6}y=#Rd7^=>&Xl3&a?wKj!;
zn~xTCZ6oq~t&=$SH}U^s$$b&l3@(($^;~{*N*m&|D~(AD*{CY4yeN?&K-zsdMZJAo
zydNi9I-kOTp0bGeyCI_AJsxfj7x<`EOugpJH$+oP^NFTbZD_4LRV>_PK$y63(ry<7
z*nT}__ed?VpBtfz_J=Z&NeTP<*GCTSObntZn2`BLIF?njxOkefQo=8s;7$!Y8%i7-
zN%_~-Nb-cZ@@&zD`cxh=&nTah-bZk^_ZGHX4jaX6tuH5+3t@CwqfnT~CAjxUj#4%q
zk9r@qP7CEg<NZ`q{v+c5Leu4E8gS7!{fnpACOXOw9GdAa&c^pF%6ZFG0yrL<zW38*
zE*fra*=`+4hr|u<5=Uac3w7ffe|-|**yO2vt5_ox*g8J3FVIHAz4#sbU$C)^dAC?p
zQ%IQFq!+(Q{EU}VIBqc(e=KI4D0|0-i~40hCh=b}&CSe(mIUAF$X`%sB=J+}!O^5n
zHl9k*FMiY(g6+Fs)w-37lxHEq6L}1H{kXIB1Bq`-k9sz!_X;2y?R`G(Ef+>V6Dsem
z(SpDFjgy>8HvAPk#5<k|pt7|=*6chNmgABZT9No&D$*(*A;E#e_tfP&eS9o&d3z#e
zuOTdgM$ewSLMMLGY8FF7AM<C9>~feOK>Z^HnfpAlJ_9v+=ebOT^@P7j2;mT>>%wUP
zO9;x^jh(+6x!9M!kKyh`$Dz{v@`7PDSfeK=pFYJ$=v3AF8|HD*V&WeYK$zI)n++BX
zt|0y)D|g(KMgeXv3AB~e=3!4>;q#{(nfPAmFV%6y01I`D-s|cT`*FAGww=X=*{MSU
zT^R;46Pt?GO0Y37ZTf|{{Up8wUtyjc&n0>KoY%Ud3@n!o)BaB4Ynf(h$|VNzzq#MU
zdp45xKfduQ)n`EJ$N3xINc<`)JE0dANaDx0@c741xLEP&VQV3mj-+>Qr_Oa|<74|o
z<I#r#tj@RDrMioY?FFBl)k_#)pHA+axS0*Bz!!VwGKHA6R{K|fJ`WO&oi3HDnBcAc
z?tN*1jnlVJjy!Y|0-IDUePy}u=%j}AE@i;4yE=BY9UGJBddA5+1^9i#Xw9-~hS+oE
z^aBMq2BLSLf5#*1>)Urb`e3{e9+nYD=EU&eMj5P{xu1!lVAGtB2RV3mbJ!?-k`M(x
zER~H1d3b6D_8nI`Hib1NR?i^gySv(9UbK+xSj4Kv8uD<6XZZ7>DiiWyxs|JJiGH}X
z-P7iqA=K;in#v%uzUMY4N)hHQ->NGvyM&DopMRP8`U=5NsxQ_O@^HRPanj;8ZM@s9
zd}0fck7Q!=_Wiwld^mEfq<R^_XZF`V++V;zV!$^kr_*e#oL0ClVU7?c)k=m`BOabO
zKM9e2Nk`Y_uxXOd*+`Ufd+I<J;Bv;&y)uiq(BGi#@j;0R`=*Zr!Uzuf-uNsuOc4;K
z+G|qXGed;REE<yaVZhs@t6%;J8_Vlr;?r*mae-|bsou<m-zS9^BU|Y39CVB5pU8pT
z;PlPqL_a&S9>2_=LDt8<>S_HKI$oC9*Mzlk;9BOM={!vUE7O&fo<u`1c+&Ux$uVKA
z{LrkvjE&0gJC`mSFNEW*2fAj&{|$#N&AnJaM{#V&^E0IXpVyb`yh|3KAmf|=qCnyw
z<33iGUSObbj&{q(U=kmGD$n`CC(pk*=e^@17cqM?Ch-YAQK^{Hl|k?$ujbaCvYA4>
zD0Uq#TFQf1?9(?}HJKQRI$PXPNc`72dUp9Q0oK0D)Qe*9z)sX^_~}i@hedfuE|B>8
zg-H?FyHNm>I?BCRR~}02S_`*LB}{6()uf;(4#sbzq^UCq{@E*?oW7Tf+l+qiRdRHM
z9&Gz!OxC|?`R>#W4P*N6<Zn;mVamroUhQ={T8_2Ue^Vg&!<I*1Cy>0f$Ncqo&vgX9
zgnsG|-#~V1SDP?gb2wOLwB_s69s%q{Ww&?>ct|o>nqouzgZ~95?XoQ#3@)y<7L^e~
z{*qFfW;qvHY8CT0oFMp3pm**wm4gLhj^Ujz_>i&_DK{Ya5Ow$R6@JpOdBLwoCSTaN
zHL}0?l(qm$zNwq;RT|=bQ7wOGFcXh|lxnJ-BL3M*B75qv0OynDq*jdZV7L-nM}9DH
zt<2JdJAs2GRUO?!*9GYI7(Z$y&qL2D!?foE4CuTpywy`q{I6@UY($_Ceae4j44M-@
zuWa&XJCZ+L-5rveew~dzMl7RcyAb-<8o#wq=VJP%-5QNybVL@pPksG?jTHk+V>LB|
zaBHY7*6}bxHEU3;Pm2ks?@2Cmv)EV?qq68gjsQC&7jLgBCjQ?f(qnBW18X#sjc(96
z_)`4J=qK?X0S)U<tQ9BnTb_1xQVRoT8%29B9AM*sMemQB2Lw>-dGR4rn~U1`ndio}
zGLYWF5qEh=>?cqE(G)EKx;6y1J^XEmotkSq(*z7`Kk%hs|28(xD#`h`{}5nAR9!u3
z5*K%4{HV3J>6r0h$(x-o*;w<SQh78+fU`f9?j?_MAuDz{zyCS|8pn;svlB=@xQkP8
zvsHlbz*|>c_i!=gmyMytQ#u&MVrlb9{JfN!bC<PNfYVk5JB|?iZrETYvS~j97yBZY
z2hZjp)<@FHYKH(Cn#oOSqC7PH*hMcR`9tXan4+RCHjHL(h}oSkK!Wy@6sLABR@U$z
z3^tHANs7~SXRIRrWmeb%zZyP%en0GVw2X^o(!sCSgpr-c+@9*a8XO1$^CV0+3vk2p
ztE#;x7lRr0#;rC4U#(5+PN^gD*_v@Evx^Vim}@HY2tNKfqlwv;&j5ZH+nP*fqksRx
z`{8v0<k`k6YLofF^+m~LzBUvR9~R1ijic%t<|xRK{3Ii)rMidU&q3|a3;9I8%wgFB
zCrJEtStB=?DS$`I=q;6}B)_Xw{(8)hj^1R~fij|>W<#Rg1qA{)B@`6Dd{5-#*I5yw
zPsi^GE2mFM<-qt};g0L%{+!Zz3vW4-{5QN(a(ONT?#d;0JN{y0c97;x-zXuTFM8}2
zy5A7(&OIUvCz17Q>b2Tj#fHLB8Ll^5h=aMtnNc>x|2kweELCKn$K|eFD!H#=`PfNf
zh~TsIB{S-a2tIsK(WOT4m9WuSF~o?2#B&e&=eiL5sqj)Wj^x*e-hNizn#91e&6{@8
zUb7KidnY0Kwg8?HD!0zx<|4GCBG&yE9q!Y5E_XU{P`Uktx|NW`$1~2R@g-dN7&X6;
ze@}<ANmvEt1{(n#^)5Dxg|IQ7-tnB^&)-Uk4|vtYer{XR-rq35)5@g%3f)9Mj+X5A
ztl(ndCF{c%!WghiSy%T}l#S!@<r9r`Nc{0jcm4d!5dGGd>wgjaR%oWXo3@g~&&tpO
zgK-4E4!@b?znljPYj=7UlZl|H@4jg}IoPE&uzD@wH_EprXdfc_Ei=u??^rYgE1P>2
z`Kt+^`=rLMVIRR)(#c<+*7H!CIVZF81QYXa%9MHPO~f6g&pCY&LL@uebJYnxl^+(L
zoKitYk79cxi>zPb?wWHUF9h(CZQXtDIKgLg!cu$f7+A4ur)wYKcO;LzmX80;hyR@w
z$LQ)@tj$zaJ{CvT-@e=VF~Mgw3JYe|zay^!^k%dtsSy7XclwQ72LpLiR94Z#$ogq4
z+oRzn1b_CToB@*mp{SxGVIc!QwXap*CHT@SJ|s`;ngEAKw#yn2UfBLY%HC6d5&Uj_
zz^i(c;O{r@-<&%xz$c?;tv8dnNEJ!WzImL11O2-Bx@<P~Pd~_=exBH`R(bVtdoC2b
zC~w$gC(DXfA#_*g!1u+X2!&t)ijMIG=e`+2+uY@ox-K34F2fr4o^mjHVW%C9BLuhW
z%dQ%dzeV#OX-#aVLq$qZvIra;$h>-C^o{`T*FKq-^ccc_^9ikoM+yJIwydnM;y_Ex
zWykFmLfp@*m$bDa`Ny0P&UHx!q}Q`sHXLT-S9;KM+iWtwW7E5ryyN1<o!4m+>j{3!
z3d>^7;ovYu&7y$JU%%crSeoGbERWdnL30?eoM*;T(<J$mQ}*L=seEYo)DD?WBJsVZ
zYH8VJ28NAILx*0mA^(ZKxwS}u#X1*0X}NL1FN{0<h}c)-?7IpOknwv+D(>Da#FXy`
zN0$=*Wfo`Q;53r|4}I8^@QL7qv~>x?w=W5xrT6wK^BE8FE($O5xeT22cAQs}!bY#n
zlh2$j0yHfU%+V(6le<}8cYzKQ=W0J!{8~%!^ZI-5&k+C12sv$E5zEDi>~CTVJ;=@?
zD{a3nmkpV-Yc^(*^_g})=W|dB$<JqbjmoU0BOrryRJw)iZ;Yf~&wngH4M$4f<FO&a
zWw;ZPp40JsNZj;NG8@tr_9J%Z2*07b@zSw#BtP}?P#Jb-pmwR}UgaPT6f5I2{7L^B
zg-Si^34ZI`7}(eC!NBW;E9cwhvGL06iAb1%04wM8p6qZX{8;|VS=w{x$hqbyv`b<m
zJx;&bfW#NK!oeB&M1Iq@4P09iz`&LFKBu?z8zA@6<%QkEe-6HxZYdBmf{uBA)?2bZ
z6Fa2T)TVGy`(5DWq$WUz?osxHI^w^+pZKmu^uKgu=R&m;1m9kdON*H)K&q^>zT#~j
z0uD!V6bG32GU(4!zt2XV<L}PgWdc0;)qAPE&JY^E6^{lkV&KJU!(G8H96YkUUm+JE
zK>G41r@C50sKg&NH#$$^+otT0<SGNaUBg?}@SfPO{?aC|LDK)Ej(NFN40Ijcdh&)Q
z8;UDlkDm4wfTEJhKCqSOha#;+VHe4-eayz6C-`8+?#~y-g_3<dea#Am)g*rD9EmC&
zV!&<v`l+p=9847%5gEQk+P5EkRFgpbujM@P7Q&BMnAl9)=uhIW>vzkwodWRZ$KM+y
z`Qgm(i>}=vJL@ZL_mw;#`fX%eYQ||5VpsoE-o<TPw9&H{r^qu<eDlZZ8cQ}TddwO>
zX9=+EvFZ1L?M8SZb$b3L!jCmA{;SWD;9oI=q#Ye({?T8bMs8wrp+^aK=qLUGg$8!M
z95xb~MP3<|ko&<`{3T9L0_{UXTIp#z&NK+0UAV80=?gRU>WF`4jMw=*aW)Uuh3nQf
zw9`qxQ=0pL#liIQk^D3%0a|w-RNrux3(6kOGIjFi@D2sZ-0(p*zQu{@xrhkSU83f<
z=^YRHUKzh!$eWJm){Q&dvx@`g+LW!%5+r|@sH(3S;3D&Ax<V$=U&Xn1@0|Wg{O9MS
z`>MqLzbdDDjEWKbIwC5{P@>~(x0>e#DT43snLhU?`>1CXDQ_2K6aS`s(YAO61CEge
z+$Tgo4TQWRS4knJUN*ftXC)66+vmxgAa7zSynLsjZpeYf8*Z?qmJrj$ouxtvzG!`U
zI{wQf2JZW2KT)O-e#=OqnL+k_TbhR-K3QXko7L9Z!vx<a9dvayJI97*)Dh*?!9t`&
z$Xko;B>aq_%Z;XX62D&iZCXI^afBzWNB)=)DhIyD?kD&fpSo8aA^1w=Li>@3sT>4x
zBy&nsgb?iU_%e_1XHKjD`XYi4c2^Aknz)#aLQ~T~KQ7_d{aq6tju?WkA~t{BO$LP1
zi!NUNW&lI$rWx}HKXHM=uN~UYMf*vGllkHdvY(K-+l}Cl7N14GeA)%rqihjnvX6(G
z@oAod5C(?U?U^P)@~;yOYg_+%EkMZ$+Zk`9dAQwsw|fVl*jK2$*j5V;5_W!Hngs!_
z81a40i2sUrIu<iR^5^hZUbX{yY`DACo5_**X%d;kR#h{?`xPn2>be*hpYlTT`fZYb
z=6#|4d1XIV^zYa4|MRaR{&SQj+w&BPIN6>T87CrX@~4=ZjD1PUfBByOKYj`D|Ew?g
z&#V8Ge||CW=jlIxT|^&F`7?yz|8@wBzlWgo_m3vL|NJV#JCghZg{aN{{Ek4J9LXTX
W|Ng%dQGEM<I{9C8>2^e=>;D29L5xKJ

literal 0
HcmV?d00001


From 04b5ccf7320eb7225e06d55877b26e0979288ca7 Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Wed, 24 Mar 2021 17:16:24 -0700
Subject: [PATCH 04/10] truncated SVD L96 example

---
 .gitignore                                    |   2 +
 examples/Lorenz/GModel.jl                     |   3 +-
 examples/Lorenz/Lorenz_example.jl             |  83 +++++++++++-------
 examples/Lorenz/Project.toml                  |   1 +
 examples/Lorenz/outputinput_output_pairs.jld2 | Bin 17704 -> 0 bytes
 src/GaussianProcessEmulator.jl                |  52 +++++++----
 src/MarkovChainMonteCarlo.jl                  |   9 +-
 src/Utilities.jl                              |  18 +++-
 8 files changed, 114 insertions(+), 54 deletions(-)
 delete mode 100644 examples/Lorenz/outputinput_output_pairs.jld2

diff --git a/.gitignore b/.gitignore
index ba418360f..83b96794c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -14,3 +14,5 @@ docs/site/
 *.jl.mem
 *.vscode*
 *.DS_Store
+*.jld2
+output*
diff --git a/examples/Lorenz/GModel.jl b/examples/Lorenz/GModel.jl
index 7bc08985a..944415f9c 100644
--- a/examples/Lorenz/GModel.jl
+++ b/examples/Lorenz/GModel.jl
@@ -5,7 +5,8 @@ using DocStringExtensions
 using Random
 using Distributions
 using LinearAlgebra
-using DifferentialEquations
+#using DifferentialEquations
+using OrdinaryDiffEq
 using FFTW
 
 export run_G
diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl
index 524bc1601..67d08dd67 100644
--- a/examples/Lorenz/Lorenz_example.jl
+++ b/examples/Lorenz/Lorenz_example.jl
@@ -23,7 +23,9 @@ using CalibrateEmulateSample.ParameterDistributionStorage
 using CalibrateEmulateSample.DataStorage
 using CalibrateEmulateSample.Observations
 
-rng_seed = 4137
+#rng_seed = 4137
+#rng_seed = 413798
+rng_seed = 2413798
 Random.seed!(rng_seed)
 
 # Output figure save directory
@@ -45,7 +47,7 @@ end
 dynamics = 2 # Transient is 2
 # Statistics integration length
 # This has to be less than 360 and 360 must be divisible by Ts_days
-Ts_days = 90. # Integration length in days
+Ts_days = 30. # Integration length in days
 # Stats type, which statistics to construct from the L96 system
 # 4 is a linear fit over a batch of length Ts_days
 # 5 is the mean over a batch of length Ts_days
@@ -91,7 +93,7 @@ end
 if dynamics == 2
 	#prior_means = [F_true+0.5, A_true+0.5]
 	prior_means = [F_true, A_true]
-	prior_stds = [2.0, 0.5*A_true]
+	prior_stds = [3.0, 0.5*A_true]
 	d1 = Parameterized(Normal(prior_means[1], prior_stds[1]))
 	d2 = Parameterized(Normal(prior_means[2], prior_stds[2]))
 	prior_distns = [d1, d2]
@@ -125,11 +127,12 @@ dt = 1/64.
 # Start of integration
 t_start = 800.
 # Data collection length
-if dynamics==1
-    T = 2.
-else
-    T = 360. / τc
-end
+#if dynamics==1
+#    T = 2.
+#else
+#    T = 360. / τc
+#end
+T = 360. / τc
 # Batch length
 Ts = 5. / τc # Nondimensionalize by L96 timescale
 # Integration length
@@ -163,6 +166,7 @@ lorenz_params = GModel.LParams(F_true, ω_true, A_true)
 gt = dropdims(GModel.run_G_ensemble(params_true, lorenz_settings), dims=2)
 
 # Compute internal variability covariance
+n_samples = 50 
 if var_prescribe==true
     n_samples = 100
     yt = zeros(length(gt),n_samples)
@@ -175,7 +179,6 @@ if var_prescribe==true
     end
 else
     println("Using truth values to compute covariance")
-    n_samples = 20 
     yt = zeros(length(gt), n_samples)
     for i in 1:n_samples
 	    lorenz_settings_local = GModel.LSettings(dynamics, stats_type, t_start+T*(i-1), 
@@ -189,12 +192,19 @@ else
     Γy = cov(yt, dims=2)
     
     println(Γy)
+    println(size(Γy))
+    println(rank(Γy))
 end
 
 
 # Construct observation object
 truth = Observations.Obs(yt, Γy, data_names)
+# Truth sample for EKP
 truth_sample = truth.mean
+#sample_ind=randperm!(collect(1:n_samples))[1]
+#truth_sample = yt[:,sample_ind]
+#println("Truth sample:")
+#println(sample_ind)
 ###
 ###  Calibrate: Ensemble Kalman Inversion
 ###
@@ -251,18 +261,24 @@ end
 ###
 
 # Emulate-sample settings
-standardize = false
-truncate_svd = 1.0
+standardize = true
+truncate_svd = 0.95
 
 gppackage = GaussianProcessEmulator.GPJL()
 pred_type = GaussianProcessEmulator.YType()
 
 # Standardize the output data
-norm_factor = get_standardizing_factors(yt) 
+# Use median over all data since all data are the same type
+yt_norm = vcat(yt...)
+norm_factor = get_standardizing_factors(yt_norm) 
 println(size(norm_factor))
-norm_factor = vcat(norm_factor...)
+#norm_factor = vcat(norm_factor...)
+norm_factor = fill(norm_factor, size(truth_sample))
+println("Standardization factors")
+println(norm_factor)
 
 # Get training points from the EKP iteration number in the second input term  
+N_iter = 5;
 input_output_pairs = Utilities.get_training_points(ekiobj, N_iter)
 normalized = true
 
@@ -276,23 +292,28 @@ gpobj = GaussianProcessEmulator.GaussianProcess(input_output_pairs, gppackage;
 
 # Check how well the Gaussian Process regression predicts on the
 # true parameters
-if log_normal==false
-    y_mean, y_var = GaussianProcessEmulator.predict(gpobj, reshape(params_true, :, 1), 
-                                   transform_to_real=true)
-else
-	y_mean, y_var = GaussianProcessEmulator.predict(gpobj, reshape(log.(params_true), :, 1), 
-                                   transform_to_real=true)
+if truncate_svd==1.0
+    if log_normal==false
+        y_mean, y_var = GaussianProcessEmulator.predict(gpobj, reshape(params_true, :, 1), 
+                                       transform_to_real=true)
+    else
+    	y_mean, y_var = GaussianProcessEmulator.predict(gpobj, reshape(log.(params_true), :, 1), 
+                                       transform_to_real=true)
+    end
+    
+    println("GP prediction on true parameters: ")
+    println(vec(y_mean))
+    println(" GP variance")
+    println(diag(y_var[1],0))
+    println("true data: ")
+    if standardize
+        println(truth.mean ./ norm_factor)
+    else
+        println(truth.mean)
+    end
+    println("GP MSE: ")
+    println(mean((truth.mean - vec(y_mean)).^2))
 end
-
-println("GP prediction on true parameters: ")
-println(vec(y_mean))
-println(" GP variance")
-println(diag(y_var[1],0))
-println("true data: ")
-println(truth.mean)
-println("GP MSE: ")
-println(mean((truth.mean - vec(y_mean)).^2))
-
 ###
 ###  Sample: Markov Chain Monte Carlo
 ###
@@ -354,11 +375,11 @@ save_directory = figure_save_directory*y_folder
 for idx in 1:n_params
     if idx == 1
         param = "F"
-	xbounds = [7.95, 8.05]
+	xbounds = [F_true-1.0, F_true+1.0]
 	xs = collect(xbounds[1]:(xbounds[2]-xbounds[1])/100:xbounds[2])
     elseif idx == 2
         param = "A"
-	xbounds = [2.45, 2.55]
+	xbounds = [A_true-1.0, A_true+1.0]
 	xs = collect(xbounds[1]:(xbounds[2]-xbounds[1])/100:xbounds[2])
     elseif idx == 3
         param = "ω"
diff --git a/examples/Lorenz/Project.toml b/examples/Lorenz/Project.toml
index c3b779a9b..0bca61625 100644
--- a/examples/Lorenz/Project.toml
+++ b/examples/Lorenz/Project.toml
@@ -6,6 +6,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
diff --git a/examples/Lorenz/outputinput_output_pairs.jld2 b/examples/Lorenz/outputinput_output_pairs.jld2
deleted file mode 100644
index 546e4eb36205cd9cabf04adee690da8059fc90e6..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 17704
zcmeIac{r5s|2I6$*thIiX2~{FW-K#<ki4vERVpRQk}cT^rL<^KDwV`1Nk!3WDN&TP
z$r7ows7TTxl%!B7&*{_Wd*9FVJMQ~;+|QrSb9~O@7}s^pbzbN4KHtmh{W`DV;^=I~
z^Ye`g@Mq1430>jKa`cV%WjRNz4)cxXvepEwj#?2B&N4AFH8Q2JIDe{`8krdhS>{5X
z-->9~ET@T-f6M=8C9uoI(OFJfQ(9C~l0p&vGmX)I3diZsb&36-qw;#&CRZ0nH&Mz?
z>QoBMNp%*v6Yd)pK%rPFicrRpD--e&As-4wk`n2=n%s*fN<<lFNxm7oZ}ImX3PskA
z^d}-ik*3fnA!C~79lb6xfI<}=YhbL54MvW9YQ`=|`P1m1^Z)!c&ybq+=X{!j>#V<P
zjP*f^$UpVRS{N%15vv1?#&mnZSTAF(|7p+k@4FO=Q+QNBm|tkXOyAHj-|)GsBbEn5
zMFm9t{mA_2h}FJ9e;#*;2#@w%5gxGmpXaFmv?Ke^_22zCher5Dn+wOjqWyE9NcZnD
zUGdLR?4R$)#4%R>Px_zhyJB^K|JWq{=P1PgH|_svl|t#2|IgXK$A2a8uLS;;z`qjs
zR|5Y^;9m*+D}jF{@UH~^mB7Ce_*Vk||C7KUo+6M+Cp-no^vBZgD3Znp_<Q@4{Q3VS
z4KYR%{D+1Zqb0_uga2I9$ocsEsjL2U=qN_+Q>8@6w__El(iH0dWHco8#`pjhWsEEM
z&*M~b{$Goob&2p#b1l^Qe|VNLR%1nYWK47vh4SCdCn6^L?^=@oZLR+!@@0&g_|GuK
z$oc<$*5c$^cFNc!#}V1eQUX^5g#Hoj7>g72KUf}7!o!RSVyry>hno4%QQX0F+n;F)
z|5s8=|A#4&R2h@nI7)@mpJvYJ{5h(dr@Z<5agkPqf1W2bMG3Jp^ryn#D=>CSRiTKO
zvTXlW0ty8bp$tVhN(3*Nqd;2DQaJ2+`G=v84uVc71ob_a#F`t1mf<RDkPPT}&s3Fx
zO0GuuBDy{lowit4^+_OW+`FG_cOA?UqpC|84%>$AoSoXnWaFi1Zs|jcG)`qH7SNjH
zAaLd6Pvpy@({`1@`?qvt8PBrpY8<ltk#LvyNo@jn5sq5512Sl@%yV+BSH<u_JqwfK
zueM@YPErk5b<iHjrAqT8P<gnddAFkuKHT2->9CqIu3wo|{@Z>6%G<Vy#%-6u<<AQy
z2A*Pp@@-OSPLmScZtP4Luw+4YU{mX)^>PSl5AvRLo(|saPi6{C8Z@1R(JzN|!FpO}
z)bl_Js}i*@YHiR4o*xYVm_2H{XQsd9^#V=2yBz2fK~=;DwWT#%9?;Ra(WN2niYy{K
zGka>MvcV`)ZrfR`h^su`qRfx_P!;+3)9$tydZx~K+Ix|Wl1h~ux6MUivngefpMy5|
zF&j!;Xc8z-`#pZ$8Wx<gzX!j*C4-zTzgnJ&Xu^a3%81e=fxGRB^&vtP1UIdJU#=$s
zeX+X7!89Eh*zGN$71NMW@y-3vD;D1PJIT#DDS`f)vrS)*Xdu3K*5|p&vPhaz;MFQK
z9!_eP69cu!Ay`HFbgVB6?nW`^3{T49n#K3M>M?Yj@*f)PpE_!LPc_Aa8l-`Y)p3(+
z?Bx-FV*1K71FSJjR}D!TuvK~JDSdsV4kn5`@F@F2L+o4i2+(!$I5ufUo{St8f6OY)
z`@zJgC^4~*6D6RzHAO_ckO_C0g~xZ@SBBBfV}+d?RZ;NH=jY5uN!)rAa8+%mHtez%
zeqp>5MfkeNEh{?|Q6s14mQ_Z<$*cC1*{LcpUugSsU9UXutrAt(w^J9QAvRNZmu29c
z6mzRIhyiiW?71%rM{RdcS<0DvQU|{WyAH;>P@#Elvx}sn4!S>X8#vo6g&jq$H#`c}
zkTY&+Nbf8~1P;fp*%YCVr0Wm+vs*<lu=lF!hc*_N`XLu>w}_#^x%Ov^yEa7byS&=|
zRRUAgGJRjBG0_rmAW%u4hV++_{mlU^oS1I6<7lD;sxDm=-*$;XuI2if$>Io<$$vUU
zk%eU^51%jh7lWJX&|g8vm>50XwRum!IINOoAIzSp4W~kd<bfV3ths%~M`<Y&1?xt~
z7Y~b~<Kx!#ho-QgpSLT-mqA4it;c@iVkSPEEYjUGQ3Ne3-^B2FEclGL4lH~qhJ_z(
zY;>M6amygJPL5o!S$O%ZQa2Oh@*R(sU8KU<Pn<SKhO`%#_HvGn7?uaR6>h!DgyvGm
zktL6ayx3E8HncF&JQVcu;BryCoWU=-smnyqPlKdO0vg8sT4NPo%S1v&8C67H6z$<9
zZ?fZ=s0fIQ`*Bko>DH2m)j2HKT-ouL<#sW&rMIcS*vi1~RVg2?v8Xt)U4wVAh=FOd
zen@0xQ<28)8F;vsff0@{FE&67p*n5b1H@RE`#h+!+Cl`2YG*5TB{N|sd-9>@87f3x
zSv&}y&BVcDi<4M6R5%`XQmJM$ku7jAbMB{Me8)HYL0uMdj%>Vpc!Y)<gS0Il?lG}F
z`mC(;VNvX<m1+I!EE72;>34EEY1pSf+}>Tm#C$*JoT}4esNJG8&OesuD?6jqO^(dh
z{iRC$UMALQUEKZ6hKi%xEtP`IS(sCDNouvR7?#Z4xPMm?6P2lZ;)?G{qs~xTvw9i}
z36F}GxStfq^5FO>Kgsy>jjPsN^%uu{ksT~kQx={V-J5avHx)m9CX@`UC*wP#wyJfG
zIN0@R#hzL${3^fdkQOC@cOO~v8YwJDpITed_*E1a4jfO&Ol08Yl%e&#>SB2HB~Le4
zlZ77>YBa)#fBhtN?m?0~3vcfy?%i}&0vZRTm<|t_I6ic9@_9EAeD>39q5G2c5-+G;
z>MDlhbeC+9^|Csj>Cvhrjx&W1?_`nwqwkHc3^0?%G8SWM{sI<;j_Svhc8jB1ZnJVT
zsV|$c(V~3?6~2o#CSS5(LH=4>aYr{5t?VhYEngD*ANiR0)slw03McJC4O!4unj~$q
zjfOK%RM+(sF)<^-G?*ixVtm}Ph+Sm8DlJCG<uy|=dDcY_K@}5PuMYMVvB~<L9r;%A
zhUg!6;}Kh;|I@w?O^ct;#ExFuik9=D7^pb4VzGdQSpf!#lJ}?>U_Ea3A7r8`Qa0$D
zkO~R+Vk-uVh13rPEGcUmuH2lUH>|;e>5;GpZ%)t<wdHDxn+OYkwe3<Gy-!8cledN<
zZA>UJkG%gSMT5*9k?dtdOsHQ8Pl;mDkd-(5em0GTNf#bp6r<1}csMzFehm|wW`?Yd
znNG#q-_jc`$$AO2nmegEH29R?G^`h8A@-NgJv|2+9#waj_7izLF>82hc$fxTZ!RO8
z=*JK1z_=GzMX~+cbQutTVQW1zP`ptTso}A+9V(ef+i;EdPKgHIG>xfKiF{qpwHa*A
z5W&OL-=a3ep2O|$E)F90)N+U0K9tSGedh5(w_BpnwSS~2p38(8mHnpu1QoaQ9FFYJ
zW5Kvx`^K?&8g^Yekvf^!yAeO{tA_y<9ZNYUS1Gfga(iLY<K;AbovPm`X2QbYyV(1B
zIy3~{Rh42i5c?64o*Oik2FdzwFU!iA=%3|SM?XfzflGl)XDno4mjElSYS6$`3v(Y-
zVnMC+M}VuM7!u_4j`))PX+yMQ&YCnNp0?>;Nv>bp{Z3b?BmT$pX^rOyk!Nvpy>Sn@
zU%Vn&sgBHt_AWbn6&1PTb-WzdEIiy1H#xqF21fhW0E!$7j2@#42Ir}GJ)f1+bb*Px
zUn;(t92Z5oOx0BuvB%Q5TAA_ZX_#Bd@M^PW!Oi9ELT?wM|But}{37;$%|=RZD4z<;
z{fsX$rA++JU!O5`7Zn*|6;meHGx5`DepOlo(Wi=iqd5~;xZLwgr+OL<PaVZ9y2yGK
z*VLNpFQGyGOp}ui(Fb?WkPMF`qE9Pth#sJ`FdTa3QK%{v&(Bjk%7}kfiPDnJ?5D!+
zxrKg_Gz%B5sQB-p(J=kSg0@FQ-~5wpCajYp`h2KKW(Ajp58`S1^bsm%WCY1iAo85z
zPC4jvl?p5WhcYA0OzfF-@#{W0DhiJ>Y-M*b;Um80@VWz{h}P`vd#l4jRdrfyjUo+4
z_ii59OZu~Te&WPJZ5qtoJEZ1lv*2g)PT<0&VWZM=n;ueMGf3SbA(00Eb!!te$@tYY
z4TnY!Qo&c<P*6hjxh03^cf5}ZSMTbdQN%v4WDWfA45T7Y;`NTJ5-hxWvN(6!Q!2K<
zIzO}b8WS-StW@Rw$b3V3ulT)ZLUw(ptidW84)dE=f%q$K<>_m+o2a-jdU4fwD;D0R
zaAF13G{{vgqs}LIVmhTS?fsxAG~`oq-;N{k#qMrfx(*dX=c@7K2AMx~GN<Yz729^s
z8_oX6#H89e8IA%`eATM|($vqy*$?v)|5`)juOQ;{>NXSU&oAd}BL0W|s=+C?k%`6$
zObHh;8mx_CGWo>+?K;@tUQH+Q!#nwb#8)Qv4K*hgD^hW4^GTUPDHdjA$?6&JqhZ1)
zr9&Y^AEFN!e(xEiBH1|o@!R`M>^gm_CzHevg(Uf@ETTVGUd`HO+CT;E<%@y0#w^^B
znz{MiS1QDsa?HPx`52!}Q5-r>!yNZZ4}6II@7SAK+*(P6(`-?WB=Jw>ABCrk-H5;S
zmS`&?{>ArV*rEp^1iv<I`7oV~&w1k4V`UvwQ0!mN%$>-B>aJh2g(WoXe{*v*iui}g
z8`jTB;}g6%yYXfb!5`BL#J$ozMUkVwp4~i$h3<K-pDs+M;dSF=`CpS+@CsSdcV?9s
z^v|Tx>@PAgUM(g65Q+b7XWq=&w~q;(`uDeVdPH${vB$F(vi_{sT}AhfP%(Ywb+Zg&
z4+m#Ae79LfLvMBI7Y!8_UWPu(7`;VBMRdg6tNBE~w<g}2J&uYsOZwi7>DTi!$+Q#1
zf87t+?3_#HSF}=HL4(M@H^SSbjKRXTis99dX3|h+6xo^IMB>8?+vVQGo{srk*%n3g
zb-u%@c|2npw#BUDn-Kl$?Oi{+_YM_QDyPVu#DlT{E93KVRIGP;(w62<@Ymf-=~LXv
z_!kXX+Y@=mDCe7`GHCGr>xZujkymC_Ps<{rpV|J+*3!g3(MxB4m{U$gVrk&UOoF#8
zZr9oi>}mKb+eh~Ru?MeIu@~bfQ&GP>zv&6FU-9i;qo0Yrz4KP`5V^xdWWtNZZ}q7t
zV|Ja1Bl@11MtiZ0^tbU*z3m&KA8kvOb&_JKSj3jsGkMFz?E42_dl7$fzB_kK7r{e|
zr|wFvCjRY}e46C~M}j|owsm%Xqv1$%w^lP*e_b(=Wy}#`9}Z>?u4MkpXWA^<K9kr}
zz__mTS4{YZbO`Sf|KTrNuRToU;ZYfFyzPM~5+~Bcl6y#hA=`H3XOZ@2{G5E{Gx47V
z&X*r_Q4t&F6DmvQ_geJ#ok9|yQ$zh@gB3`;o#wnR<~bFSx9fEV$a)>N`29R}R1~{w
z91ispdHRYO)P~j5FiGmor8=Syj;@}<U#c{0a;wz&LgL5l&eO_Ui9g+ulPD!i^sCn)
zM(Y>B6ExaZr?Zn-Fkmgxx3;0dvm;ZoM4W}KDI9xx1P$C=t;_CeER?L?^SY;p3Ljtj
z9s#*N9ap@WdY20AiY1IPQvXR!(SrTTG@Ogi(@rAund@P{{p&d@l6*Ck$J?+lUi7+F
z(P<ir4rd*0A^xbGKka}j!3Q0|n+Bhe_!+EuuQ=$kD2BW!(_8o~WTbvt*K9yT+~sk(
zgTy~wa>xi#jTD3FLZ=4bR3-+BzgI+?(eO;<PWWDOo%|{}cA7dBzc0k-`+gwyQLJ{S
zllZ^m%=uHh$owRw92bX?_HI@8x~`Wc{xoo*^w2gEZ({Tw3J5+`klWXEmrupJ3C~7m
zlKMXB{j3LzY3MTZqW{VCH!;OT|K?ZzOQt_Y!i<&wZKm(YB$+;#UgE}28unV@ch~h`
z3wJ&wO@2jdT+l&_l=z_mCJ&Xd!mxXnb+NGj>=}G#VnX)LmP(f)+s)qR>1}pOsMxS-
zse88#mK;jCcyp&YMo(PI8(cF9dN0<jUbaL)H2?QtjlLOXs9fx^C^toyLrYs~AqQ!P
zdv1$eF+xS@+9$cofrI;RJX-R^5;|JR=YsAS;c%7Q;8RmGd@*iX;^n4`3RREj*JqR=
zo4Yv(AFkSR^|lS1N~XZDf1|aFoC%tr8bv0Tagch4vf7hqV&Hpa4^Eyw{HLGHwJV;C
zC;8q}%eI(6zNPaqOHUIE9@c(YpvFPGy}nDrO&zcrq`pTFYT=i@T2t)-B^)p{mt8$?
zJeJO-{|-HFi7LOCIWw&Jc-Vcsxa+nqT1$h2?0b!{DNQ$;KV*aLwm#iepZGA8xjyjv
zg&O85ERZYwZipN&-by!BONiy%J2_fwf(iP$36(0^crv|4+&yt3;+>a9&`w&xI&K}W
z%Ebf{kteQsY#)a~dy$&zvj*5-wq$l$h6R>+8CZD@7-6sGxp2=6RcPG0G_2vn#(O*6
zp;hzw=nQMQA!j%VcT=-e{cns%NXDay_lI<$cOcDX^Xn-nop7`NkeM0YeTvN1Gc`cg
z<{$iOhAD(28+KW2GzM+z$4ZY@HmLHJCT*`Y5IQS4b>eGXgngQ~+ON+P9M!$%CBAIT
zu~d|fqlsYdm3dq1-ZG%bxZD`gAVim8+v>#M6EXkr`Ce(EA|`$|nU%Ck7rXD@`U{2*
zxGE{e^KiC?dv{2u<zWu2Q@>D7NSmQ7Th?r}dMc!f4Ej7LPsWA#r86V+c#v`BUQb^q
z#1Y4PnwmC%_L>ZrtXymC4j5MyZOg}UiLd-jM+*#`jZf%2WQv2aHS<l6a`5WXaJPrD
zGU~Q1_c)|W7%9$B)09aT$hi`BNBW~N<X_jc?ex>Yvy^RvIocfPu6m#DwA2JCGfRU<
zH|ZmR=hy4FO%_w;zxO)8A{q3$E8<n3^zq7?>!xL?jWm01d(VR(w!t&acbl7y!;H3k
z?;Exj7-`cyzMpG~y6yMNEsy8|hh)CUL~(Kc!l%ppr{<t1?u<>+Gr`=e>u!E)q+{)2
z=E0-4Igr<iWUm#N!0VNa+Og^Skg>5%$R-x9{nO#5Uw}4bU5|aZvc(kdwc-p@z4Wl`
zgWL*J{Ing-uJn%bQ$e9ySo^6Y3wX7&57kr{!%jroCUmws-kHXB_3Cq=T<)73aLEKW
z0w#p+PvhX&cI&Eac{)0+${fv|^fCO+p3}XA4_|@A%x^S(Y_1DlxMi&<n$;Jbdg84H
z{%+CshQF*a`*X^P{O{&?(P(@1);Iy^8PhjybTvnArD}TpPcr`-Ki(^8n!;{>U5J*p
zKC1hdoLV>A2#1E}Jum3D1}~*!Jd-0pf6kI?qsbGnx&6v-9~~1Ey^c$s8Y9G@J}0oQ
zj|=}9O+g#dHBqizQqWgv0RM(57cvH{AhXRe?9z231O?=_ZXJ=t<c9&>x)yp+f4)#{
zgQ6a$N9}k^-^oBnLES*&kr7)z%QI6QDi!g3@s_TG`%K`ru$4yH#D*-FopbG*I$;{5
z+P2ToMbE*oyYp%^usQzL18rwT^qpVQ=u|pn%TQ~7)>tcp3p1x&T<ocfxf;frcOGj(
zMCXu#NLhz%g4HE;&vpg$df;O@rodNMxysQ_fOYFDl#b}?!tCw+hqF6OuoH)TALy9E
zu}pA(rTHYt{-}|fTc?dnGi~H@nS@cxZS=YwLF6~n{M-ozKK52kew;Gb04u+l6!D^s
zG3(_<F&CT32xK)}@mXn#6se8Zn)6ubizwd^!s226(~!pc%acfktFY#gmkw@ENK(7l
z(``FNL)9Adl#poX#*cetjT6c7M_)yl;n>ulM=TS0IC7IN`BhVZ$wk(#Lu<?+@jGp;
z@i=3g>SN~D+mXyQXS#1#7tznr0{8rIKAcr5=bW9!#$8_U*69t3@HE`^mSd;`hGz*)
zC4|iT=epf@E*fI+qv`sRdODghdwlykHm3Bgj{f-73{zgjhfVy&gW9;!m!{U*IFVs`
zO2V82*{oZIdm;spVb!&niyIOq<oZtP1Pw%<-CrakYJh~(<8#{|2r*n?a3uUME=*db
zn^Uh5JAAt~T_lnXTK2pbWgE?K<4Id~l?fM5)Cz|0Zr4PytZ~U;o&nCT?GkS~z=yR~
zjNifcY=jyG*_TTwU}oH4=xd5LG^opM0t)!(8_=yVn#%=iqP%{{GdePjlJ?DeZ-Dy5
zxsK6fep{2@PTn|~hjy!zb@Tfe*z>N*=gvtEPDv%6@!ZA7$Cf1p#cvI<dSKvK<$MO_
z#s`+0br~S!>6?l$iWzd&)b}f8^04!{rcN-Aj%h}}oQK`mu-AAgdE8Nm76bW~7pJ+<
z+5bv958BwC6R)=A1{+7>l*D%&6Cmt?yv?})E;e|(&VMk%K>gYE3%u5{;rv5AY@eDC
zPv1Q|uql!U*-_OLK`R}yMx~kYhuBysQu^_-g#fcNJFm@h;Nn_;rp(C(I{5QiZ|jZY
zz~#NyBG!FASSmEOvNspsUGFYmHk*lw&T<9W)*SR2q#GR95un-dw^pgBAto0}-h8L3
z4Tp{O1uN4H(6ZPlp)!GwW`~><>sN-TyXH4Rlpv?V(p6~zyV=-94J;p^2~pJfS!73y
zA+(&oAGc4{#J0&FBX1NN;K@d+=PwH({CB>Z9~I0)tE0G%+^9D4CA4B{j5+Z8^6A``
zWIl|y?uhnqnuI)jI504*g)1vB7lpR!!*t2I=;fbHabxwz)8}#~Az;US`-731ct3ft
zCY#(pQFv$j_89_LhDkhg51#~)mlZ!gdTHSChSf@@jrs_4;U-2J3$Tc~H%(28hg8Wq
zN972!L`ybIY<tVb-n1<TO<MTKyEa~7e!n3errg<?&Zpy<p3?06tp-?k(Pq)&cmcFc
zeg{0+z=iw+ho43RS~zom?eX7wY_uo2?Y*NXM5xr6y1_0(Y-fAzlf12ki476O3X*Kp
zJa;iY|4;zcNkRKo5WJNttyF$gnt`Vxb%)FyIM5Wc3b=18!0i>EoTqzm5pNXU=9s~N
z<X-1Y-&w@I^0Vj9n=8PzvuhGqt%j&yCoSDqr3Iy^C6WEd*tk)ZwO9PE0G(g5<#*<A
zF~QC1tA)8XuKVx0@x|W&4tIicntTLUl9aV1R-B7Z<t+J>csfouh_71Z&c-O=M+3U~
zI8DhKzilNK+NqD5?9%D5l?pKQQsLl-jl0qcfe_nUTfRSR;=*b8?E06{+BoSvS9Y2y
z8}Ih-Ug0yt48y*Ci#Jtskxo+(y>d+rdhb`y&sJiiHhg?Tw2>KVFC0i|Si^(2+Y6bb
zx!Q1ysv24<!@>6z&dx$pA(Vzf0z%ph34`S#F@d9rqiKag8zR4Q-&d0#OADc`N6BiV
z^3YhZXvccO+*~>n9-H%lgZ`C_tR?<JXc+59bi{F?aYKAk{bx<wYk%Et@6AS}dFi3u
zvSxte)*oWmd00_X>Kdh^h4P!HFW%hG#uJB$@8XXN@U&P%p(fK1*LKwyb@wn3{Dt}9
zXek?O>t{aEc*V!y(3=g>v6C=??lXVxat#DLaKKH1C$tZ__;WhVVb<I6P9&Fyvd>{I
z!&Pd~U2eon4l;m_^r@P?d1QQ-8pk*K@}T;APo!d?CPIdMUcbJ_#>Vq<PZtF8v1ys>
z*UK?n^!uIr{Wgy<3;Iqr!Ui@T?>es6Li9T+`^ZauV=jJ_*{$UVXu)k=pw%U^zV)Xa
zca=>sLw~<jz3l=cNLfeLMh$Dj`BioId{qu=uADidK`{qq$4_fhbz?Lt+}vHBJsz7c
zOp>YeV<XdMhUrMP5bLgfk#MQtVdT`kzKT!U*l8ix%pv~Uf2Qm4kTU`_evsfRxAI_8
zRM)?1923zx7uP$U<=}*MV%oT2Av}%7e{)pf!G1;giZ#Y$r=~w*lhXzc=5Mx{tC1!^
zrBwsda6S)??Z%Z^%jnqKr=5{>g$?etJq_W8W-w&fkM_UkVr%V|M}j;CTGY>Myh_&h
zhpa)~xduLVE>~vlHsyl5>IqXuO&hzloDBAp=hts3UV6}y=#Rd7^=>&Xl3&a?wKj!;
zn~xTCZ6oq~t&=$SH}U^s$$b&l3@(($^;~{*N*m&|D~(AD*{CY4yeN?&K-zsdMZJAo
zydNi9I-kOTp0bGeyCI_AJsxfj7x<`EOugpJH$+oP^NFTbZD_4LRV>_PK$y63(ry<7
z*nT}__ed?VpBtfz_J=Z&NeTP<*GCTSObntZn2`BLIF?njxOkefQo=8s;7$!Y8%i7-
zN%_~-Nb-cZ@@&zD`cxh=&nTah-bZk^_ZGHX4jaX6tuH5+3t@CwqfnT~CAjxUj#4%q
zk9r@qP7CEg<NZ`q{v+c5Leu4E8gS7!{fnpACOXOw9GdAa&c^pF%6ZFG0yrL<zW38*
zE*fra*=`+4hr|u<5=Uac3w7ffe|-|**yO2vt5_ox*g8J3FVIHAz4#sbU$C)^dAC?p
zQ%IQFq!+(Q{EU}VIBqc(e=KI4D0|0-i~40hCh=b}&CSe(mIUAF$X`%sB=J+}!O^5n
zHl9k*FMiY(g6+Fs)w-37lxHEq6L}1H{kXIB1Bq`-k9sz!_X;2y?R`G(Ef+>V6Dsem
z(SpDFjgy>8HvAPk#5<k|pt7|=*6chNmgABZT9No&D$*(*A;E#e_tfP&eS9o&d3z#e
zuOTdgM$ewSLMMLGY8FF7AM<C9>~feOK>Z^HnfpAlJ_9v+=ebOT^@P7j2;mT>>%wUP
zO9;x^jh(+6x!9M!kKyh`$Dz{v@`7PDSfeK=pFYJ$=v3AF8|HD*V&WeYK$zI)n++BX
zt|0y)D|g(KMgeXv3AB~e=3!4>;q#{(nfPAmFV%6y01I`D-s|cT`*FAGww=X=*{MSU
zT^R;46Pt?GO0Y37ZTf|{{Up8wUtyjc&n0>KoY%Ud3@n!o)BaB4Ynf(h$|VNzzq#MU
zdp45xKfduQ)n`EJ$N3xINc<`)JE0dANaDx0@c741xLEP&VQV3mj-+>Qr_Oa|<74|o
z<I#r#tj@RDrMioY?FFBl)k_#)pHA+axS0*Bz!!VwGKHA6R{K|fJ`WO&oi3HDnBcAc
z?tN*1jnlVJjy!Y|0-IDUePy}u=%j}AE@i;4yE=BY9UGJBddA5+1^9i#Xw9-~hS+oE
z^aBMq2BLSLf5#*1>)Urb`e3{e9+nYD=EU&eMj5P{xu1!lVAGtB2RV3mbJ!?-k`M(x
zER~H1d3b6D_8nI`Hib1NR?i^gySv(9UbK+xSj4Kv8uD<6XZZ7>DiiWyxs|JJiGH}X
z-P7iqA=K;in#v%uzUMY4N)hHQ->NGvyM&DopMRP8`U=5NsxQ_O@^HRPanj;8ZM@s9
zd}0fck7Q!=_Wiwld^mEfq<R^_XZF`V++V;zV!$^kr_*e#oL0ClVU7?c)k=m`BOabO
zKM9e2Nk`Y_uxXOd*+`Ufd+I<J;Bv;&y)uiq(BGi#@j;0R`=*Zr!Uzuf-uNsuOc4;K
z+G|qXGed;REE<yaVZhs@t6%;J8_Vlr;?r*mae-|bsou<m-zS9^BU|Y39CVB5pU8pT
z;PlPqL_a&S9>2_=LDt8<>S_HKI$oC9*Mzlk;9BOM={!vUE7O&fo<u`1c+&Ux$uVKA
z{LrkvjE&0gJC`mSFNEW*2fAj&{|$#N&AnJaM{#V&^E0IXpVyb`yh|3KAmf|=qCnyw
z<33iGUSObbj&{q(U=kmGD$n`CC(pk*=e^@17cqM?Ch-YAQK^{Hl|k?$ujbaCvYA4>
zD0Uq#TFQf1?9(?}HJKQRI$PXPNc`72dUp9Q0oK0D)Qe*9z)sX^_~}i@hedfuE|B>8
zg-H?FyHNm>I?BCRR~}02S_`*LB}{6()uf;(4#sbzq^UCq{@E*?oW7Tf+l+qiRdRHM
z9&Gz!OxC|?`R>#W4P*N6<Zn;mVamroUhQ={T8_2Ue^Vg&!<I*1Cy>0f$Ncqo&vgX9
zgnsG|-#~V1SDP?gb2wOLwB_s69s%q{Ww&?>ct|o>nqouzgZ~95?XoQ#3@)y<7L^e~
z{*qFfW;qvHY8CT0oFMp3pm**wm4gLhj^Ujz_>i&_DK{Ya5Ow$R6@JpOdBLwoCSTaN
zHL}0?l(qm$zNwq;RT|=bQ7wOGFcXh|lxnJ-BL3M*B75qv0OynDq*jdZV7L-nM}9DH
zt<2JdJAs2GRUO?!*9GYI7(Z$y&qL2D!?foE4CuTpywy`q{I6@UY($_Ceae4j44M-@
zuWa&XJCZ+L-5rveew~dzMl7RcyAb-<8o#wq=VJP%-5QNybVL@pPksG?jTHk+V>LB|
zaBHY7*6}bxHEU3;Pm2ks?@2Cmv)EV?qq68gjsQC&7jLgBCjQ?f(qnBW18X#sjc(96
z_)`4J=qK?X0S)U<tQ9BnTb_1xQVRoT8%29B9AM*sMemQB2Lw>-dGR4rn~U1`ndio}
zGLYWF5qEh=>?cqE(G)EKx;6y1J^XEmotkSq(*z7`Kk%hs|28(xD#`h`{}5nAR9!u3
z5*K%4{HV3J>6r0h$(x-o*;w<SQh78+fU`f9?j?_MAuDz{zyCS|8pn;svlB=@xQkP8
zvsHlbz*|>c_i!=gmyMytQ#u&MVrlb9{JfN!bC<PNfYVk5JB|?iZrETYvS~j97yBZY
z2hZjp)<@FHYKH(Cn#oOSqC7PH*hMcR`9tXan4+RCHjHL(h}oSkK!Wy@6sLABR@U$z
z3^tHANs7~SXRIRrWmeb%zZyP%en0GVw2X^o(!sCSgpr-c+@9*a8XO1$^CV0+3vk2p
ztE#;x7lRr0#;rC4U#(5+PN^gD*_v@Evx^Vim}@HY2tNKfqlwv;&j5ZH+nP*fqksRx
z`{8v0<k`k6YLofF^+m~LzBUvR9~R1ijic%t<|xRK{3Ii)rMidU&q3|a3;9I8%wgFB
zCrJEtStB=?DS$`I=q;6}B)_Xw{(8)hj^1R~fij|>W<#Rg1qA{)B@`6Dd{5-#*I5yw
zPsi^GE2mFM<-qt};g0L%{+!Zz3vW4-{5QN(a(ONT?#d;0JN{y0c97;x-zXuTFM8}2
zy5A7(&OIUvCz17Q>b2Tj#fHLB8Ll^5h=aMtnNc>x|2kweELCKn$K|eFD!H#=`PfNf
zh~TsIB{S-a2tIsK(WOT4m9WuSF~o?2#B&e&=eiL5sqj)Wj^x*e-hNizn#91e&6{@8
zUb7KidnY0Kwg8?HD!0zx<|4GCBG&yE9q!Y5E_XU{P`Uktx|NW`$1~2R@g-dN7&X6;
ze@}<ANmvEt1{(n#^)5Dxg|IQ7-tnB^&)-Uk4|vtYer{XR-rq35)5@g%3f)9Mj+X5A
ztl(ndCF{c%!WghiSy%T}l#S!@<r9r`Nc{0jcm4d!5dGGd>wgjaR%oWXo3@g~&&tpO
zgK-4E4!@b?znljPYj=7UlZl|H@4jg}IoPE&uzD@wH_EprXdfc_Ei=u??^rYgE1P>2
z`Kt+^`=rLMVIRR)(#c<+*7H!CIVZF81QYXa%9MHPO~f6g&pCY&LL@uebJYnxl^+(L
zoKitYk79cxi>zPb?wWHUF9h(CZQXtDIKgLg!cu$f7+A4ur)wYKcO;LzmX80;hyR@w
z$LQ)@tj$zaJ{CvT-@e=VF~Mgw3JYe|zay^!^k%dtsSy7XclwQ72LpLiR94Z#$ogq4
z+oRzn1b_CToB@*mp{SxGVIc!QwXap*CHT@SJ|s`;ngEAKw#yn2UfBLY%HC6d5&Uj_
zz^i(c;O{r@-<&%xz$c?;tv8dnNEJ!WzImL11O2-Bx@<P~Pd~_=exBH`R(bVtdoC2b
zC~w$gC(DXfA#_*g!1u+X2!&t)ijMIG=e`+2+uY@ox-K34F2fr4o^mjHVW%C9BLuhW
z%dQ%dzeV#OX-#aVLq$qZvIra;$h>-C^o{`T*FKq-^ccc_^9ikoM+yJIwydnM;y_Ex
zWykFmLfp@*m$bDa`Ny0P&UHx!q}Q`sHXLT-S9;KM+iWtwW7E5ryyN1<o!4m+>j{3!
z3d>^7;ovYu&7y$JU%%crSeoGbERWdnL30?eoM*;T(<J$mQ}*L=seEYo)DD?WBJsVZ
zYH8VJ28NAILx*0mA^(ZKxwS}u#X1*0X}NL1FN{0<h}c)-?7IpOknwv+D(>Da#FXy`
zN0$=*Wfo`Q;53r|4}I8^@QL7qv~>x?w=W5xrT6wK^BE8FE($O5xeT22cAQs}!bY#n
zlh2$j0yHfU%+V(6le<}8cYzKQ=W0J!{8~%!^ZI-5&k+C12sv$E5zEDi>~CTVJ;=@?
zD{a3nmkpV-Yc^(*^_g})=W|dB$<JqbjmoU0BOrryRJw)iZ;Yf~&wngH4M$4f<FO&a
zWw;ZPp40JsNZj;NG8@tr_9J%Z2*07b@zSw#BtP}?P#Jb-pmwR}UgaPT6f5I2{7L^B
zg-Si^34ZI`7}(eC!NBW;E9cwhvGL06iAb1%04wM8p6qZX{8;|VS=w{x$hqbyv`b<m
zJx;&bfW#NK!oeB&M1Iq@4P09iz`&LFKBu?z8zA@6<%QkEe-6HxZYdBmf{uBA)?2bZ
z6Fa2T)TVGy`(5DWq$WUz?osxHI^w^+pZKmu^uKgu=R&m;1m9kdON*H)K&q^>zT#~j
z0uD!V6bG32GU(4!zt2XV<L}PgWdc0;)qAPE&JY^E6^{lkV&KJU!(G8H96YkUUm+JE
zK>G41r@C50sKg&NH#$$^+otT0<SGNaUBg?}@SfPO{?aC|LDK)Ej(NFN40Ijcdh&)Q
z8;UDlkDm4wfTEJhKCqSOha#;+VHe4-eayz6C-`8+?#~y-g_3<dea#Am)g*rD9EmC&
zV!&<v`l+p=9847%5gEQk+P5EkRFgpbujM@P7Q&BMnAl9)=uhIW>vzkwodWRZ$KM+y
z`Qgm(i>}=vJL@ZL_mw;#`fX%eYQ||5VpsoE-o<TPw9&H{r^qu<eDlZZ8cQ}TddwO>
zX9=+EvFZ1L?M8SZb$b3L!jCmA{;SWD;9oI=q#Ye({?T8bMs8wrp+^aK=qLUGg$8!M
z95xb~MP3<|ko&<`{3T9L0_{UXTIp#z&NK+0UAV80=?gRU>WF`4jMw=*aW)Uuh3nQf
zw9`qxQ=0pL#liIQk^D3%0a|w-RNrux3(6kOGIjFi@D2sZ-0(p*zQu{@xrhkSU83f<
z=^YRHUKzh!$eWJm){Q&dvx@`g+LW!%5+r|@sH(3S;3D&Ax<V$=U&Xn1@0|Wg{O9MS
z`>MqLzbdDDjEWKbIwC5{P@>~(x0>e#DT43snLhU?`>1CXDQ_2K6aS`s(YAO61CEge
z+$Tgo4TQWRS4knJUN*ftXC)66+vmxgAa7zSynLsjZpeYf8*Z?qmJrj$ouxtvzG!`U
zI{wQf2JZW2KT)O-e#=OqnL+k_TbhR-K3QXko7L9Z!vx<a9dvayJI97*)Dh*?!9t`&
z$Xko;B>aq_%Z;XX62D&iZCXI^afBzWNB)=)DhIyD?kD&fpSo8aA^1w=Li>@3sT>4x
zBy&nsgb?iU_%e_1XHKjD`XYi4c2^Aknz)#aLQ~T~KQ7_d{aq6tju?WkA~t{BO$LP1
zi!NUNW&lI$rWx}HKXHM=uN~UYMf*vGllkHdvY(K-+l}Cl7N14GeA)%rqihjnvX6(G
z@oAod5C(?U?U^P)@~;yOYg_+%EkMZ$+Zk`9dAQwsw|fVl*jK2$*j5V;5_W!Hngs!_
z81a40i2sUrIu<iR^5^hZUbX{yY`DACo5_**X%d;kR#h{?`xPn2>be*hpYlTT`fZYb
z=6#|4d1XIV^zYa4|MRaR{&SQj+w&BPIN6>T87CrX@~4=ZjD1PUfBByOKYj`D|Ew?g
z&#V8Ge||CW=jlIxT|^&F`7?yz|8@wBzlWgo_m3vL|NJV#JCghZg{aN{{Ek4J9LXTX
W|Ng%dQGEM<I{9C8>2^e=>;D29L5xKJ

diff --git a/src/GaussianProcessEmulator.jl b/src/GaussianProcessEmulator.jl
index 49026e6c2..3f39115b0 100644
--- a/src/GaussianProcessEmulator.jl
+++ b/src/GaussianProcessEmulator.jl
@@ -21,7 +21,7 @@ export predict
 
 export GPJL, SKLJL
 export YType, FType
-export svd_transform, svd_reverse_transform_mean_cov, get_standardizing_factors
+export svd_transform, svd_reverse_transform_mean_cov
 
 """
     GaussianProcessesPackage
@@ -106,15 +106,18 @@ function GaussianProcess(
 
     # Initialize vector of GP models
     models = Any[]
-    # Number of models (We are fitting one model per output dimension)
-    N_models = output_dim
    
 
     # TO DO: Standardize the data here
     # Can use the full time median or some user define function?
     if standardize
         #norm_factors = get_standarizing_factors() 
+	println(size(get_outputs(input_output_pairs)))
+	output_values = get_outputs(input_output_pairs) ./ norm_factor
 	obs_noise_cov = obs_noise_cov ./ (norm_factor .* norm_factor')
+	println(size(output_values))
+    else
+	output_values = get_outputs(input_output_pairs)
     end
 
     # Normalize the inputs if normalized==true 
@@ -123,16 +126,31 @@ function GaussianProcess(
     if normalized
         # Normalize (NB the inputs have to be of) size [input_dim × N_samples] to pass to GPE())
         sqrt_inv_input_cov = sqrt(inv(Symmetric(cov(get_inputs(input_output_pairs), dims=2))))
-        GPinputs = sqrt_inv_input_cov * (get_inputs(input_output_pairs).-input_mean) 
+        GPinputs = sqrt_inv_input_cov * (get_inputs(input_output_pairs).-input_mean)
     else
         GPinputs = get_inputs(input_output_pairs)
     end
     
     
-    # Transform data if obs_noise_cov available (if obs_noise_cov==nothing, transformed_data is equal to data)
-    transformed_data, decomposition = svd_transform(get_outputs(input_output_pairs), 
+    # Transform data if obs_noise_cov available 
+    # (if obs_noise_cov==nothing, transformed_data is equal to data)
+    transformed_data, decomposition = svd_transform(output_values, 
 						    obs_noise_cov, truncate_svd=truncate_svd)
+    println(size(transformed_data))
+    #println(transformed_data)
+
+    # Overwrite the input-output pairs because of the size discrepency
+    if truncate_svd<1.0
+	# Overwrite the input_output_pairs structure
+	input_output_pairs_old = deepcopy(input_output_pairs)
+	input_output_pairs = PairedDataContainer(get_inputs(input_output_pairs_old),
+					      transformed_data, data_are_columns=true)
+        input_dim, output_dim = size(input_output_pairs, 1)
+    end
 
+    # Number of models (We are fitting one model per output dimension)
+    N_models = output_dim; #size(transformed_data)[1]
+    
     # Use a default kernel unless a kernel was supplied to GaussianProcess
     if GPkernel==nothing
         println("Using default squared exponential kernel, learning length scale and variance parameters")
@@ -244,7 +262,10 @@ function GaussianProcess(
     # Can use the full time median or some user define function?
     if standardize
         #norm_factors = get_standarizing_factors() 
+	output_values = get_outputs(input_output_pairs) ./ norm_factor
 	obs_noise_cov = obs_noise_cov ./ (norm_factor .* norm_factor')
+    else
+	output_values = get_outputs(input_output_pairs)
     end
 
     
@@ -263,7 +284,7 @@ function GaussianProcess(
 
     # Transform data if obs_noise_cov available (if obs_noise_cov==nothing, 
     # transformed_data is equal to data)
-    transformed_data, decomposition = svd_transform(get_outputs(input_output_pairs), 
+    transformed_data, decomposition = svd_transform(output_values, 
 						    obs_noise_cov, truncate_svd=truncate_svd)
 
     if GPkernel==nothing
@@ -504,10 +525,10 @@ function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, No
             # Apply truncated SVD
             n = size(obs_noise_cov)[1]
             decomposition.S[k+1:n] = zeros(n-k)
-            sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(SVD_local.S))
+            sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S))
             transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
             transformed_data = transformed_data[1:k, :];
-            transformed_data = convert(Matrix{Float64},transformed_data');
+            #transformed_data = convert(Matrix{Float64},transformed_data');
 	else
             decomposition = svd(obs_noise_cov)
             sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
@@ -539,10 +560,10 @@ function svd_transform(data::Vector{FT},
             # Apply truncated SVD
             n = size(obs_noise_cov)[1]
             decomposition.S[k+1:n] = zeros(n-k)
-            sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(SVD_local.S))
+            sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S))
             transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
-            transformed_data = transformed_data[1:k, :];
-            transformed_data = convert(Matrix{Float64},transformed_data');
+            transformed_data = transformed_data[1:k];
+	    #transformed_data = permutedims(transformed_data, [2, 1]);
 	else
             decomposition = svd(obs_noise_cov)
             sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
@@ -588,12 +609,5 @@ function svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::Array{FT, 2}, dec
     return transformed_μ, transformed_σ2
 end
 
-function get_standardizing_factors(data::Array{FT,2}) where {FT}
-    # Input: data size: N_data x N_ensembles
-    # Ensemble median of the data
-    norm_factor = median(data,dims=2)
-    return norm_factor
-end
-
 
 end 
diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl
index 56b2794d0..fe35ee5e0 100644
--- a/src/MarkovChainMonteCarlo.jl
+++ b/src/MarkovChainMonteCarlo.jl
@@ -75,7 +75,7 @@ function MCMC(
     max_iter::IT,
     algtype::String,
     burnin::IT,
-    norm_factors::Union{Array{FT, 1}, Nothing};
+    norm_factor::Union{Array{FT, 1}, Nothing};
     svdflag=true,
     standardize=false,
     truncate_svd=1.0) where {FT<:AbstractFloat, IT<:Int}
@@ -84,11 +84,15 @@ function MCMC(
     param_init_copy = deepcopy(param_init)
     
     # Standardize MCMC input?
+    println(obs_sample)
+    println(obs_noise_cov)
     if standardize
-        obs_sample = obs_sample ./ norm_factors
+        obs_sample = obs_sample ./ norm_factor;
 	cov_norm_factor = norm_factor .* norm_factor;
 	obs_noise_cov = obs_noise_cov ./ cov_norm_factor;
     end
+    println(obs_sample)
+    println(obs_noise_cov)
 
     # We need to transform obs_sample into the correct space 
     if svdflag
@@ -97,6 +101,7 @@ function MCMC(
     else
         println("Assuming independent outputs.")
     end
+    println(obs_sample)
     
     # first row is param_init
     posterior = zeros(length(param_init_copy),max_iter + 1)
diff --git a/src/Utilities.jl b/src/Utilities.jl
index c657a8562..24dbf4248 100644
--- a/src/Utilities.jl
+++ b/src/Utilities.jl
@@ -12,7 +12,7 @@ export get_training_points
 export get_obs_sample
 export orig2zscore
 export zscore2orig
-
+export get_standardizing_factors
 
 """
     extract_training_points(ekp::EnsembleKalmanProcess{FT, IT, P}, N_train_iterations::IT) where {FT,IT, P}
@@ -113,4 +113,20 @@ function zscore2orig(Z::AbstractMatrix{FT},
     return X
 end
 
+function get_standardizing_factors(data::Array{FT,2}) where {FT}
+    # Input: data size: N_data x N_ensembles
+    # Ensemble median of the data
+    norm_factor = median(data,dims=2)
+    return norm_factor
+end
+
+function get_standardizing_factors(data::Array{FT,1}) where {FT}
+    # Input: data size: N_data*N_ensembles (splatted)
+    # Ensemble median of the data
+    norm_factor = median(data)
+    return norm_factor
+end
+
+
+
 end # module

From e8bde68fcdbd58f2e80f3f3bcd5abbb5d3704755 Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Tue, 30 Mar 2021 16:40:08 -0700
Subject: [PATCH 05/10] removed DiffEq dependency

---
 examples/Lorenz/Project.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/examples/Lorenz/Project.toml b/examples/Lorenz/Project.toml
index 0bca61625..db6900f9d 100644
--- a/examples/Lorenz/Project.toml
+++ b/examples/Lorenz/Project.toml
@@ -1,6 +1,5 @@
 [deps]
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"

From 494638d42ea4f837306cade644879f39563c8af4 Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Wed, 31 Mar 2021 10:59:16 -0700
Subject: [PATCH 06/10] Fixed diagonal matrix issue in log_gp

---
 examples/Lorenz/GModel.jl         |  2 +-
 examples/Lorenz/Lorenz_example.jl | 16 ++++++++++++++++
 src/GaussianProcessEmulator.jl    | 31 ++++++++++++++++++-------------
 src/MarkovChainMonteCarlo.jl      | 27 ++++++++++++++++-----------
 src/Utilities.jl                  | 16 ----------------
 5 files changed, 51 insertions(+), 41 deletions(-)

diff --git a/examples/Lorenz/GModel.jl b/examples/Lorenz/GModel.jl
index 944415f9c..c7b5d4491 100644
--- a/examples/Lorenz/GModel.jl
+++ b/examples/Lorenz/GModel.jl
@@ -6,7 +6,7 @@ using Random
 using Distributions
 using LinearAlgebra
 #using DifferentialEquations
-using OrdinaryDiffEq
+#using OrdinaryDiffEq
 using FFTW
 
 export run_G
diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl
index 67d08dd67..bf9b9db91 100644
--- a/examples/Lorenz/Lorenz_example.jl
+++ b/examples/Lorenz/Lorenz_example.jl
@@ -28,6 +28,22 @@ using CalibrateEmulateSample.Observations
 rng_seed = 2413798
 Random.seed!(rng_seed)
 
+
+function get_standardizing_factors(data::Array{FT,2}) where {FT}
+    # Input: data size: N_data x N_ensembles
+    # Ensemble median of the data
+    norm_factor = median(data,dims=2)
+    return norm_factor
+end
+
+function get_standardizing_factors(data::Array{FT,1}) where {FT}
+    # Input: data size: N_data*N_ensembles (splatted)
+    # Ensemble median of the data
+    norm_factor = median(data)
+    return norm_factor
+end
+
+
 # Output figure save directory
 example_directory = @__DIR__
 println(example_directory)
diff --git a/src/GaussianProcessEmulator.jl b/src/GaussianProcessEmulator.jl
index 3f39115b0..deed88d88 100644
--- a/src/GaussianProcessEmulator.jl
+++ b/src/GaussianProcessEmulator.jl
@@ -591,19 +591,24 @@ elements on the main diagonal (i.e., the variances), we return the full
 covariance at each point, as a vector of length N_predicted_points, where 
 each element is a matrix of size output_dim × output_dim
 """
-function svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::Array{FT, 2}, decomposition::SVD) where {FT}
-
-    output_dim, N_predicted_points = size(σ2)
-    # We created meanvGP = D_inv * Vt * mean_v so meanv = V * D * meanvGP
-    sqrt_singular_values= Diagonal(sqrt.(decomposition.S))
-    transformed_μ = decomposition.V * sqrt_singular_values * μ
-
-    transformed_σ2 = [zeros(output_dim, output_dim) for i in 1:N_predicted_points]
-    # Back transformation
-
-    for j in 1:N_predicted_points
-        σ2_j = decomposition.V * sqrt_singular_values * Diagonal(σ2[:,j]) * sqrt_singular_values * decomposition.Vt
-        transformed_σ2[j] = σ2_j
+function svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::Array{FT, 2}, 
+                                        decomposition::SVD;
+					truncate_svd::FT=1.0) where {FT}
+    if truncate_svd==1.0
+        output_dim, N_predicted_points = size(σ2)
+        # We created meanvGP = D_inv * Vt * mean_v so meanv = V * D * meanvGP
+        sqrt_singular_values= Diagonal(sqrt.(decomposition.S))
+        transformed_μ = decomposition.V * sqrt_singular_values * μ
+
+        transformed_σ2 = [zeros(output_dim, output_dim) for i in 1:N_predicted_points]
+        # Back transformation
+
+        for j in 1:N_predicted_points
+            σ2_j = decomposition.V * sqrt_singular_values * Diagonal(σ2[:,j]) * sqrt_singular_values * decomposition.Vt
+            transformed_σ2[j] = σ2_j
+        end
+    else
+        println("To do")
     end
 
     return transformed_μ, transformed_σ2
diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl
index fe35ee5e0..d8018eede 100644
--- a/src/MarkovChainMonteCarlo.jl
+++ b/src/MarkovChainMonteCarlo.jl
@@ -150,9 +150,10 @@ function get_posterior(mcmc::MCMC)
     
 end
 
-function mcmc_sample!(mcmc::MCMC{FT}, g::Vector{FT}, gvar::Vector{FT}) where {FT}
+function mcmc_sample!(mcmc::MCMC{FT}, g::Vector{FT}, 
+                      gcov::Union{Matrix{FT},Diagonal{FT}}) where {FT}
     if mcmc.algtype == "rwm"
-        log_posterior = log_likelihood(mcmc, g, gvar) + log_prior(mcmc)
+        log_posterior = log_likelihood(mcmc, g, gcov) + log_prior(mcmc)
     end
 
     if mcmc.log_posterior[1] isa Nothing # do an accept step.
@@ -175,6 +176,10 @@ function mcmc_sample!(mcmc::MCMC{FT}, g::Vector{FT}, gvar::Vector{FT}) where {FT
 
 end
 
+function mcmc_sample!(mcmc::MCMC{FT}, g::Vector{FT}, gvar::Vector{FT}) where {FT}
+    return mcmc_sample!(mcmc, g, Diagonal(gvar))
+end
+
 function accept_ratio(mcmc::MCMC{FT}) where {FT}
     return convert(FT, mcmc.accept[1]) / mcmc.iter[1]
 end
@@ -182,23 +187,23 @@ end
 
 function log_likelihood(mcmc::MCMC{FT},
                         g::Vector{FT},
-                        gvar::Vector{FT}) where {FT}
+                        gcov::Union{Matrix{FT},Diagonal{FT}}) where {FT}
     log_rho = FT[0]
-    if gvar == nothing
-        diff = g - mcmc.obs_sample
-        log_rho[1] = -FT(0.5) * diff' * (mcmc.obs_noise_cov \ diff)
-    else
+    #if gcov == nothing
+    #    diff = g - mcmc.obs_sample
+    #    log_rho[1] = -FT(0.5) * diff' * (mcmc.obs_noise_cov \ diff)
+    #else
 	# det(log(Γ))
 	# Ill-posed numerically for ill-conditioned covariance matrices with det≈0
         #log_gpfidelity = -FT(0.5) * log(det(Diagonal(gvar))) # = -0.5 * sum(log.(gvar))
 	# Well-posed numerically for ill-conditioned covariance matrices with det≈0
-	full_cov = Diagonal(gvar)
-	eigs = eigvals(full_cov)
+	#full_cov = Diagonal(gvar)
+	eigs = eigvals(gcov)
 	log_gpfidelity = -FT(0.5) * sum(log.(eigs))
 	# Combine got log_rho
 	diff = g - mcmc.obs_sample
-        log_rho[1] = -FT(0.5) * diff' * (Diagonal(gvar) \ diff) + log_gpfidelity
-    end
+        log_rho[1] = -FT(0.5) * diff' * (gcov \ diff) + log_gpfidelity
+    #end
     return log_rho[1]
 end
 
diff --git a/src/Utilities.jl b/src/Utilities.jl
index 24dbf4248..da2dfe662 100644
--- a/src/Utilities.jl
+++ b/src/Utilities.jl
@@ -12,7 +12,6 @@ export get_training_points
 export get_obs_sample
 export orig2zscore
 export zscore2orig
-export get_standardizing_factors
 
 """
     extract_training_points(ekp::EnsembleKalmanProcess{FT, IT, P}, N_train_iterations::IT) where {FT,IT, P}
@@ -113,20 +112,5 @@ function zscore2orig(Z::AbstractMatrix{FT},
     return X
 end
 
-function get_standardizing_factors(data::Array{FT,2}) where {FT}
-    # Input: data size: N_data x N_ensembles
-    # Ensemble median of the data
-    norm_factor = median(data,dims=2)
-    return norm_factor
-end
-
-function get_standardizing_factors(data::Array{FT,1}) where {FT}
-    # Input: data size: N_data*N_ensembles (splatted)
-    # Ensemble median of the data
-    norm_factor = median(data)
-    return norm_factor
-end
-
-
 
 end # module

From ef2935c4678a07bd9dcf2b5856cf1d092265a958 Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Wed, 31 Mar 2021 11:01:57 -0700
Subject: [PATCH 07/10] removed diffeq

---
 examples/Lorenz/Manifest.toml | 339 ----------------------------------
 examples/Lorenz/Project.toml  |   1 -
 2 files changed, 340 deletions(-)

diff --git a/examples/Lorenz/Manifest.toml b/examples/Lorenz/Manifest.toml
index 575e711eb..73976a39e 100644
--- a/examples/Lorenz/Manifest.toml
+++ b/examples/Lorenz/Manifest.toml
@@ -1,34 +1,17 @@
 # This file is machine-generated - editing it directly is not advised
 
-[[AbstractAlgebra]]
-deps = ["InteractiveUtils", "LinearAlgebra", "Markdown", "Random", "RandomExtensions", "SparseArrays", "Test"]
-git-tree-sha1 = "b773e53c62456c2937fef7e1ddd81a6bb421d81a"
-uuid = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
-version = "0.14.0"
-
 [[AbstractFFTs]]
 deps = ["LinearAlgebra"]
 git-tree-sha1 = "051c95d6836228d120f5f4b984dd5aba1624f716"
 uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
 version = "0.5.0"
 
-[[AbstractTrees]]
-git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
-uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
-version = "0.3.4"
-
 [[Adapt]]
 deps = ["LinearAlgebra"]
 git-tree-sha1 = "ffcfa2d345aaee0ef3d8346a073d5dd03c983ebe"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 version = "3.2.0"
 
-[[ArnoldiMethod]]
-deps = ["LinearAlgebra", "Random", "StaticArrays"]
-git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
-uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
-version = "0.1.0"
-
 [[Arpack]]
 deps = ["Arpack_jll", "Libdl", "LinearAlgebra"]
 git-tree-sha1 = "2ff92b71ba1747c5fdd541f8fc87736d82f40ec9"
@@ -65,32 +48,15 @@ git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9"
 uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
 version = "1.0.0"
 
-[[BandedMatrices]]
-deps = ["ArrayLayouts", "Compat", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"]
-git-tree-sha1 = "8c83ee44dc9835263760ad4e77ed4eed4b3490c1"
-uuid = "aae01518-5342-5314-be14-df237901396f"
-version = "0.15.25"
-
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
-[[BoundaryValueDiffEq]]
-deps = ["BandedMatrices", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "NLsolve", "Reexport", "SparseArrays"]
-git-tree-sha1 = "fe34902ac0c3a35d016617ab7032742865756d7d"
-uuid = "764a87c0-6b3e-53db-9096-fe964310641d"
-version = "2.7.1"
-
 [[Bzip2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "c3598e525718abcc440f69cc6d5f60dda0a1b61e"
 uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
 version = "1.0.6+5"
 
-[[CEnum]]
-git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
-uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
-version = "0.4.1"
-
 [[Cairo_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
 git-tree-sha1 = "e2f47f6d8337369411569fd45ae5753ca10394c6"
@@ -145,16 +111,6 @@ git-tree-sha1 = "ac5f2213e56ed8a34a3dd2f681f4df1166b34929"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
 version = "0.12.6"
 
-[[Combinatorics]]
-git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
-uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
-version = "1.0.2"
-
-[[CommonSolve]]
-git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f"
-uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
-version = "0.2.0"
-
 [[CommonSubexpressions]]
 deps = ["MacroTools", "Test"]
 git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
@@ -179,12 +135,6 @@ git-tree-sha1 = "6231e40619c15148bcb80aa19d731e629877d762"
 uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 version = "1.5.1"
 
-[[ConstructionBase]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "48920211c95a6da1914a06c44ec94be70e84ffff"
-uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-version = "1.1.0"
-
 [[Contour]]
 deps = ["StaticArrays"]
 git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
@@ -234,52 +184,10 @@ version = "0.4.13"
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
-[[DelayDiffEq]]
-deps = ["DataStructures", "DiffEqBase", "LinearAlgebra", "Logging", "NonlinearSolve", "OrdinaryDiffEq", "Printf", "RecursiveArrayTools", "Reexport", "UnPack"]
-git-tree-sha1 = "0c7f367792397f754728a3db960d7ef1eb3d8f09"
-uuid = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
-version = "5.29.1"
-
 [[DelimitedFiles]]
 deps = ["Mmap"]
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 
-[[DiffEqBase]]
-deps = ["ArrayInterface", "ChainRulesCore", "DataStructures", "DocStringExtensions", "FunctionWrappers", "IterativeSolvers", "LabelledArrays", "LinearAlgebra", "Logging", "MuladdMacro", "NonlinearSolve", "Parameters", "Printf", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "ZygoteRules"]
-git-tree-sha1 = "f5d290a67a2105b4cf531b0e3b5178106a661e39"
-uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.57.9"
-
-[[DiffEqCallbacks]]
-deps = ["DataStructures", "DiffEqBase", "ForwardDiff", "LinearAlgebra", "NLsolve", "OrdinaryDiffEq", "RecipesBase", "RecursiveArrayTools", "StaticArrays"]
-git-tree-sha1 = "d4c4a3f442ab749b6b895c514b0be984c75d6d67"
-uuid = "459566f4-90b8-5000-8ac3-15dfb0a30def"
-version = "2.16.0"
-
-[[DiffEqFinancial]]
-deps = ["DiffEqBase", "DiffEqNoiseProcess", "LinearAlgebra", "Markdown", "RandomNumbers"]
-git-tree-sha1 = "db08e0def560f204167c58fd0637298e13f58f73"
-uuid = "5a0ffddc-d203-54b0-88ba-2c03c0fc2e67"
-version = "2.4.0"
-
-[[DiffEqJump]]
-deps = ["ArrayInterface", "Compat", "DataStructures", "DiffEqBase", "FunctionWrappers", "LinearAlgebra", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "StaticArrays", "TreeViews", "UnPack"]
-git-tree-sha1 = "3ec8d5eeb792897b28ef79a851d834ce1415498f"
-uuid = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
-version = "6.13.0"
-
-[[DiffEqNoiseProcess]]
-deps = ["DiffEqBase", "Distributions", "LinearAlgebra", "Optim", "PoissonRandom", "QuadGK", "Random", "Random123", "RandomNumbers", "RecipesBase", "RecursiveArrayTools", "Requires", "ResettableStacks", "StaticArrays", "Statistics"]
-git-tree-sha1 = "f300e85c99c79bdd0434df097d2803ae36e6ccca"
-uuid = "77a26b50-5914-5dd7-bc55-306e6241c503"
-version = "5.6.0"
-
-[[DiffEqPhysics]]
-deps = ["DiffEqBase", "DiffEqCallbacks", "ForwardDiff", "LinearAlgebra", "Printf", "Random", "RecipesBase", "RecursiveArrayTools", "Reexport", "StaticArrays"]
-git-tree-sha1 = "8f23c6f36f6a6eb2cbd6950e28ec7c4b99d0e4c9"
-uuid = "055956cb-9e8b-5191-98cc-73ae4a59e68a"
-version = "3.9.0"
-
 [[DiffResults]]
 deps = ["StaticArrays"]
 git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
@@ -292,18 +200,6 @@ git-tree-sha1 = "214c3fcac57755cfda163d91c58893a8723f93e9"
 uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
 version = "1.0.2"
 
-[[DifferentialEquations]]
-deps = ["BoundaryValueDiffEq", "DelayDiffEq", "DiffEqBase", "DiffEqCallbacks", "DiffEqFinancial", "DiffEqJump", "DiffEqNoiseProcess", "DiffEqPhysics", "DimensionalPlotRecipes", "LinearAlgebra", "MultiScaleArrays", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "RecursiveArrayTools", "Reexport", "SteadyStateDiffEq", "StochasticDiffEq", "Sundials"]
-git-tree-sha1 = "221b9a427fc8970be5b65171c25f0a6ba0e1f394"
-uuid = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
-version = "6.16.0"
-
-[[DimensionalPlotRecipes]]
-deps = ["LinearAlgebra", "RecipesBase"]
-git-tree-sha1 = "af883a26bbe6e3f5f778cb4e1b81578b534c32a6"
-uuid = "c619ae07-58cd-5f6d-b883-8f17bd6a98f9"
-version = "1.2.0"
-
 [[Distances]]
 deps = ["LinearAlgebra", "Statistics"]
 git-tree-sha1 = "366715149014943abd71aa647a07a43314158b2d"
@@ -350,17 +246,6 @@ git-tree-sha1 = "1402e52fcda25064f51c77a9655ce8680b76acf0"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
 version = "2.2.7+6"
 
-[[ExponentialUtilities]]
-deps = ["LinearAlgebra", "Printf", "Requires", "SparseArrays"]
-git-tree-sha1 = "712cb5af8db62836913970ee035a5fa742986f00"
-uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18"
-version = "1.8.1"
-
-[[ExprTools]]
-git-tree-sha1 = "10407a39b87f29d47ebaca8edbc75d7c302ff93e"
-uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
-version = "0.1.3"
-
 [[FFMPEG]]
 deps = ["FFMPEG_jll", "x264_jll"]
 git-tree-sha1 = "9a73ffdc375be61b0e4516d83d880b265366fe1f"
@@ -385,11 +270,6 @@ git-tree-sha1 = "5a0d4b6a22a34d17d53543bd124f4b08ed78e8b0"
 uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
 version = "3.3.9+7"
 
-[[FastClosures]]
-git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef"
-uuid = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
-version = "0.3.2"
-
 [[FastGaussQuadrature]]
 deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"]
 git-tree-sha1 = "5829b25887e53fb6730a9df2ff89ed24baa6abf6"
@@ -444,11 +324,6 @@ git-tree-sha1 = "0d20aed5b14dd4c9a2453c1b601d08e1149679cc"
 uuid = "559328eb-81f9-559d-9380-de523a88c83c"
 version = "1.0.5+6"
 
-[[FunctionWrappers]]
-git-tree-sha1 = "241552bc2209f0fa068b6415b1942cc0aa486bcc"
-uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
-version = "1.1.2"
-
 [[Future]]
 deps = ["Random"]
 uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
@@ -517,11 +392,6 @@ git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
 uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 version = "0.1.0"
 
-[[Inflate]]
-git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
-uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
-version = "0.1.2"
-
 [[IniFile]]
 deps = ["Test"]
 git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
@@ -555,12 +425,6 @@ git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
 uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 version = "1.3.0"
 
-[[IterativeSolvers]]
-deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
-git-tree-sha1 = "6f5ef3206d9dc6510a8b8e2334b96454a2ade590"
-uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
-version = "0.9.0"
-
 [[IteratorInterfaceExtensions]]
 git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
 uuid = "82899510-4779-5014-852e-03e436cf321d"
@@ -612,12 +476,6 @@ git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
 uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 version = "1.2.1"
 
-[[LabelledArrays]]
-deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"]
-git-tree-sha1 = "5e288800819c323de5897fa6d5a002bdad54baf7"
-uuid = "2ee39098-c373-598a-b85f-a56591580800"
-version = "1.5.0"
-
 [[Latexify]]
 deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
 git-tree-sha1 = "fbc08b5a78e264ba3d19da90b36ce1789ca67a40"
@@ -685,12 +543,6 @@ git-tree-sha1 = "f879ae9edbaa2c74c922e8b85bb83cc84ea1450b"
 uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
 version = "2.34.0+7"
 
-[[LightGraphs]]
-deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
-git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
-uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
-version = "1.3.5"
-
 [[LineSearches]]
 deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
 git-tree-sha1 = "f27132e551e959b3667d8c93eae90973225032dd"
@@ -710,12 +562,6 @@ git-tree-sha1 = "3242a8f411e19eda9adc49d0b877681975c11375"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
 version = "0.8.26"
 
-[[METIS_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "2dc1a9fc87e57e32b1fc186db78811157b30c118"
-uuid = "d00139f3-1899-568f-a2f0-47f597d42d70"
-version = "5.1.0+5"
-
 [[MKL_jll]]
 deps = ["IntelOpenMP_jll", "Libdl", "Pkg"]
 git-tree-sha1 = "eb540ede3aabb8284cb482aa41d00d6ca850b1f8"
@@ -758,23 +604,6 @@ version = "0.4.5"
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
-[[ModelingToolkit]]
-deps = ["ArrayInterface", "ConstructionBase", "DataStructures", "DiffEqBase", "DiffEqJump", "DiffRules", "Distributed", "Distributions", "DocStringExtensions", "IfElse", "LabelledArrays", "Latexify", "Libdl", "LightGraphs", "LinearAlgebra", "MacroTools", "NaNMath", "NonlinearSolve", "RecursiveArrayTools", "Reexport", "Requires", "RuntimeGeneratedFunctions", "SafeTestsets", "SciMLBase", "Serialization", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicUtils", "Symbolics", "TreeViews", "UnPack", "Unitful"]
-git-tree-sha1 = "5354faf8d8799a5fdc5a72fba4817b64ade8b8c0"
-uuid = "961ee093-0014-501f-94e3-6117800e7a78"
-version = "5.13.6"
-
-[[MuladdMacro]]
-git-tree-sha1 = "c6190f9a7fc5d9d5915ab29f2134421b12d24a68"
-uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
-version = "0.2.2"
-
-[[MultiScaleArrays]]
-deps = ["DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "OrdinaryDiffEq", "Random", "RecursiveArrayTools", "SparseDiffTools", "Statistics", "StochasticDiffEq", "TreeViews"]
-git-tree-sha1 = "258f3be6770fe77be8870727ba9803e236c685b8"
-uuid = "f9640e96-87f6-5992-9c3b-0743c6a49ffa"
-version = "1.8.1"
-
 [[MultivariateStats]]
 deps = ["Arpack", "LinearAlgebra", "SparseArrays", "Statistics", "StatsBase"]
 git-tree-sha1 = "8d958ff1854b166003238fe191ec34b9d592860a"
@@ -787,12 +616,6 @@ git-tree-sha1 = "50608f411a1e178e0129eab4110bd56efd08816f"
 uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
 version = "7.8.0"
 
-[[NLsolve]]
-deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"]
-git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1"
-uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
-version = "4.5.1"
-
 [[NNlib]]
 deps = ["ChainRulesCore", "Compat", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
 git-tree-sha1 = "ab1d43fead2ecb9aa5ae460d3d547c2cf8d89461"
@@ -810,12 +633,6 @@ git-tree-sha1 = "9afd724797039125e8e2cc362098f01dab60bc3a"
 uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 version = "0.4.8"
 
-[[NonlinearSolve]]
-deps = ["ArrayInterface", "FiniteDiff", "ForwardDiff", "IterativeSolvers", "LinearAlgebra", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "UnPack"]
-git-tree-sha1 = "ef18e47df4f3917af35be5e5d7f5d97e8a83b0ec"
-uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-version = "0.3.8"
-
 [[Observables]]
 git-tree-sha1 = "3469ef96607a6b9a1e89e54e6f23401073ed3126"
 uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
@@ -868,12 +685,6 @@ git-tree-sha1 = "4fa2ba51070ec13fcc7517db714445b4ab986bdf"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 version = "1.4.0"
 
-[[OrdinaryDiffEq]]
-deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "ExponentialUtilities", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "NLsolve", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
-git-tree-sha1 = "d22a75b8ae5b77543c4e1f8eae1ff01ce1f64453"
-uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-version = "5.52.2"
-
 [[PCRE_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "1b556ad51dceefdbf30e86ffa8f528b73c7df2bb"
@@ -886,12 +697,6 @@ git-tree-sha1 = "95a4038d1011dfdbde7cecd2ad0ac411e53ab1bc"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
 version = "0.10.1"
 
-[[ParameterizedFunctions]]
-deps = ["DataStructures", "DiffEqBase", "Latexify", "LinearAlgebra", "ModelingToolkit", "Reexport", "SciMLBase"]
-git-tree-sha1 = "3610913402be3856074668741326d82d02cbba5a"
-uuid = "65888b18-ceab-5e60-b2b9-181511a3b968"
-version = "5.9.0"
-
 [[Parameters]]
 deps = ["OrderedCollections", "UnPack"]
 git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345"
@@ -932,12 +737,6 @@ git-tree-sha1 = "40031283dbb5ae602aa132f423daedfc18c82069"
 uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 version = "1.11.0"
 
-[[PoissonRandom]]
-deps = ["Random", "Statistics", "Test"]
-git-tree-sha1 = "44d018211a56626288b5d3f8c6497d28c26dc850"
-uuid = "e409e4f3-bfea-5376-8464-e040bb5c01ab"
-version = "0.4.0"
-
 [[PooledArrays]]
 deps = ["DataAPI", "Future"]
 git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
@@ -992,24 +791,6 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
-[[Random123]]
-deps = ["Libdl", "Random", "RandomNumbers"]
-git-tree-sha1 = "7c6710c8198fd4444b5eb6a3840b7d47bd3593c5"
-uuid = "74087812-796a-5b5d-8853-05524746bad3"
-version = "1.3.1"
-
-[[RandomExtensions]]
-deps = ["Random", "SparseArrays"]
-git-tree-sha1 = "062986376ce6d394b23d5d90f01d81426113a3c9"
-uuid = "fb686558-2515-59ef-acaa-46db3789a887"
-version = "0.4.3"
-
-[[RandomNumbers]]
-deps = ["Random", "Requires"]
-git-tree-sha1 = "441e6fc35597524ada7f85e13df1f4e10137d16f"
-uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
-version = "1.4.0"
-
 [[Ratios]]
 git-tree-sha1 = "37d210f612d70f3f7d57d488cb3b6eff56ad4e41"
 uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
@@ -1026,18 +807,6 @@ git-tree-sha1 = "c4d54a78e287de7ec73bbc928ce5eb3c60f80b24"
 uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
 version = "0.3.1"
 
-[[RecursiveArrayTools]]
-deps = ["ArrayInterface", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
-git-tree-sha1 = "271a36e18c8806332b7bd0f57e50fcff0d428b11"
-uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "2.11.0"
-
-[[RecursiveFactorization]]
-deps = ["LinearAlgebra", "LoopVectorization", "VectorizationBase"]
-git-tree-sha1 = "a135a7d12b7b1743b0d0564a16e5e8e0585d87e5"
-uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
-version = "0.1.8"
-
 [[Reexport]]
 git-tree-sha1 = "57d8440b0c7d98fc4f889e478e80f268d534c9d5"
 uuid = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -1049,12 +818,6 @@ git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
 version = "1.1.3"
 
-[[ResettableStacks]]
-deps = ["StaticArrays"]
-git-tree-sha1 = "622b3e491fb0a85fbfeed6f17dc320a9f46d8929"
-uuid = "ae5879a3-cd67-5da8-be7f-38c6eb64a37b"
-version = "1.1.0"
-
 [[Rmath]]
 deps = ["Random", "Rmath_jll"]
 git-tree-sha1 = "86c5647b565873641538d8f812c04e4c9dbeb370"
@@ -1067,12 +830,6 @@ git-tree-sha1 = "1b7bf41258f6c5c9c31df8c1ba34c1fc88674957"
 uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
 version = "0.2.2+2"
 
-[[RuntimeGeneratedFunctions]]
-deps = ["ExprTools", "SHA", "Serialization"]
-git-tree-sha1 = "e02f14dfe3a8d3b8fc92ca80c1882bfdbc015e07"
-uuid = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
-version = "0.5.1"
-
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 
@@ -1088,18 +845,6 @@ git-tree-sha1 = "67ae90a18aa8c22bf159318300e765fbd89ddf6e"
 uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
 version = "0.5.5"
 
-[[SafeTestsets]]
-deps = ["Test"]
-git-tree-sha1 = "36ebc5622c82eb9324005cc75e7e2cc51181d181"
-uuid = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-version = "0.0.1"
-
-[[SciMLBase]]
-deps = ["ArrayInterface", "CommonSolve", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"]
-git-tree-sha1 = "9dd06ee4a1238265ddf512851682cb3db326a8a8"
-uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "1.8.5"
-
 [[ScikitLearn]]
 deps = ["Compat", "Conda", "DataFrames", "Distributed", "IterTools", "LinearAlgebra", "MacroTools", "Parameters", "Printf", "PyCall", "Random", "ScikitLearnBase", "SparseArrays", "StatsBase", "VersionParsing"]
 git-tree-sha1 = "bf4c558a0d0f1129b23e80f1ca772a9db457df4d"
@@ -1127,12 +872,6 @@ version = "1.2.16"
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
-[[Setfield]]
-deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
-git-tree-sha1 = "d5640fc570fb1b6c54512f0bd3853866bd298b3e"
-uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
-version = "0.7.0"
-
 [[SharedArrays]]
 deps = ["Distributed", "Mmap", "Random", "Serialization"]
 uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
@@ -1143,12 +882,6 @@ git-tree-sha1 = "ee010d8f103468309b8afac4abb9be2e18ff1182"
 uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
 version = "0.3.2"
 
-[[SimpleTraits]]
-deps = ["InteractiveUtils", "MacroTools"]
-git-tree-sha1 = "daf7aec3fe3acb2131388f93a4c409b8c7f62226"
-uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
-version = "0.9.3"
-
 [[Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
@@ -1162,12 +895,6 @@ version = "0.3.1"
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
-[[SparseDiffTools]]
-deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "VertexSafeGraphs"]
-git-tree-sha1 = "d05bc362e3fa1b0e2361594a706fc63ffbd140f3"
-uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
-version = "1.13.0"
-
 [[SpecialFunctions]]
 deps = ["ChainRulesCore", "OpenSpecFun_jll"]
 git-tree-sha1 = "5919936c0e92cff40e57d0ddf0ceb667d42e5902"
@@ -1202,18 +929,6 @@ git-tree-sha1 = "7c88f73b759ae16a4a1af8bda7a9563471a7baf5"
 uuid = "f3b207a7-027a-5e70-b257-86293d7955fd"
 version = "0.14.19"
 
-[[SteadyStateDiffEq]]
-deps = ["DiffEqBase", "DiffEqCallbacks", "LinearAlgebra", "NLsolve", "Reexport", "SciMLBase"]
-git-tree-sha1 = "9b908c7a5933b8ec7c2da44a477f74127baa2ce9"
-uuid = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
-version = "1.6.1"
-
-[[StochasticDiffEq]]
-deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqJump", "DiffEqNoiseProcess", "FillArrays", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Logging", "MuladdMacro", "NLsolve", "OrdinaryDiffEq", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
-git-tree-sha1 = "3109ce733c907b941eea4345b0644308e2c6da2d"
-uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
-version = "6.33.1"
-
 [[StructArrays]]
 deps = ["Adapt", "DataAPI", "Tables"]
 git-tree-sha1 = "26ea43b4be7e919a2390c3c0f824e7eb4fc19a0a"
@@ -1230,36 +945,6 @@ version = "1.5.0"
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
 uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
-[[SuiteSparse_jll]]
-deps = ["Libdl", "METIS_jll", "OpenBLAS_jll", "Pkg"]
-git-tree-sha1 = "4a2295b63d67e6f13a0b539c935ccbf218fa1143"
-uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
-version = "5.4.0+9"
-
-[[Sundials]]
-deps = ["CEnum", "DataStructures", "DiffEqBase", "Libdl", "LinearAlgebra", "Logging", "Reexport", "SparseArrays", "Sundials_jll"]
-git-tree-sha1 = "28db9e1a8fdd1b8e95cee064a6c2066897cf39c5"
-uuid = "c3572dad-4567-51f8-b174-8c6c989267f4"
-version = "4.4.1"
-
-[[Sundials_jll]]
-deps = ["CompilerSupportLibraries_jll", "Libdl", "OpenBLAS_jll", "Pkg", "SuiteSparse_jll"]
-git-tree-sha1 = "013ff4504fc1d475aa80c63b455b6b3a58767db2"
-uuid = "fb77eaff-e24c-56d4-86b1-d163f2edb164"
-version = "5.2.0+1"
-
-[[SymbolicUtils]]
-deps = ["AbstractAlgebra", "AbstractTrees", "Combinatorics", "ConstructionBase", "DataStructures", "IfElse", "LabelledArrays", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "TimerOutputs"]
-git-tree-sha1 = "5b4e14c896c9cab17807388fd2a280780ac5f258"
-uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
-version = "0.9.2"
-
-[[Symbolics]]
-deps = ["AbstractAlgebra", "DiffRules", "Distributions", "DocStringExtensions", "IfElse", "Latexify", "Libdl", "LinearAlgebra", "MacroTools", "NaNMath", "RecipesBase", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "SymbolicUtils", "TreeViews"]
-git-tree-sha1 = "b21d35c61b10de6c07753792ee721e527cba4109"
-uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-version = "0.1.10"
-
 [[TableOperations]]
 deps = ["SentinelArrays", "Tables", "Test"]
 git-tree-sha1 = "a7cf690d0ac3f5b53dd09b5d613540b230233647"
@@ -1282,24 +967,12 @@ version = "1.4.1"
 deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[[TimerOutputs]]
-deps = ["Printf"]
-git-tree-sha1 = "32cdbe6cd2d214c25a0b88f985c9e0092877c236"
-uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
-version = "0.5.8"
-
 [[TranscodingStreams]]
 deps = ["Random", "Test"]
 git-tree-sha1 = "7c53c35547de1c5b9d46a4797cf6d8253807108c"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
 version = "0.9.5"
 
-[[TreeViews]]
-deps = ["Test"]
-git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6"
-uuid = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7"
-version = "0.3.0"
-
 [[UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
@@ -1312,12 +985,6 @@ version = "1.0.2"
 [[Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
-[[Unitful]]
-deps = ["ConstructionBase", "LinearAlgebra", "Random"]
-git-tree-sha1 = "fdfbea79b5b9a305bf226eb4730321f603281290"
-uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
-version = "1.6.0"
-
 [[VectorizationBase]]
 deps = ["CpuId", "Libdl", "LinearAlgebra"]
 git-tree-sha1 = "03e2fbb479a1ea350398195b6fbf439bae0f8260"
@@ -1329,12 +996,6 @@ git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
 uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
 version = "1.2.0"
 
-[[VertexSafeGraphs]]
-deps = ["LightGraphs"]
-git-tree-sha1 = "b9b450c99a3ca1cc1c6836f560d8d887bcbe356e"
-uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f"
-version = "0.1.2"
-
 [[Wayland_jll]]
 deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
 git-tree-sha1 = "dc643a9b774da1c2781413fd7b6dcd2c56bb8056"
diff --git a/examples/Lorenz/Project.toml b/examples/Lorenz/Project.toml
index db6900f9d..f79d7ed84 100644
--- a/examples/Lorenz/Project.toml
+++ b/examples/Lorenz/Project.toml
@@ -5,7 +5,6 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"

From 59434805dcb6c9e9dc19bbc53a0d3e72fed1d9e2 Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Thu, 1 Apr 2021 15:47:30 -0700
Subject: [PATCH 08/10] added inverse SVD transform

---
 examples/Lorenz/Lorenz_example.jl |  4 ++--
 src/GaussianProcessEmulator.jl    | 35 ++++++++++++++++++-------------
 2 files changed, 22 insertions(+), 17 deletions(-)

diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl
index bf9b9db91..9404b26da 100644
--- a/examples/Lorenz/Lorenz_example.jl
+++ b/examples/Lorenz/Lorenz_example.jl
@@ -308,7 +308,7 @@ gpobj = GaussianProcessEmulator.GaussianProcess(input_output_pairs, gppackage;
 
 # Check how well the Gaussian Process regression predicts on the
 # true parameters
-if truncate_svd==1.0
+#if truncate_svd==1.0
     if log_normal==false
         y_mean, y_var = GaussianProcessEmulator.predict(gpobj, reshape(params_true, :, 1), 
                                        transform_to_real=true)
@@ -329,7 +329,7 @@ if truncate_svd==1.0
     end
     println("GP MSE: ")
     println(mean((truth.mean - vec(y_mean)).^2))
-end
+#end
 ###
 ###  Sample: Markov Chain Monte Carlo
 ###
diff --git a/src/GaussianProcessEmulator.jl b/src/GaussianProcessEmulator.jl
index deed88d88..ce9d8545a 100644
--- a/src/GaussianProcessEmulator.jl
+++ b/src/GaussianProcessEmulator.jl
@@ -46,6 +46,16 @@ abstract type PredictionType end
 struct YType <: PredictionType end
 struct FType <: PredictionType end
 
+
+# SVD decompsotion structure
+struct decomp_struct{FT<:AbstractFloat, IT<:Int}
+    V::Array{FT,2}
+    Vt::Array{FT,2}
+    S::Array{FT}
+    N::IT
+end
+
+
 """
     GaussianProcess{FT<:AbstractFloat}
 
@@ -66,7 +76,7 @@ struct GaussianProcess{FT<:AbstractFloat, GPM}
     "the Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs"
     models::Vector
     "the singular value decomposition of obs_noise_cov, such that obs_noise_cov = decomposition.U * Diagonal(decomposition.S) * decomposition.Vt."
-    decomposition::Union{SVD, Nothing}
+    decomposition::Union{SVD, decomp_struct, Nothing}
     "whether to fit GP models on normalized inputs ((inputs - input_mean) * sqrt_inv_input_cov)"
     normalized::Bool
     "prediction type (`y` to predict the data, `f` to predict the latent function)"
@@ -492,7 +502,6 @@ function predict(gp::GaussianProcess{FT, SKLJL}, new_inputs::Array{FT, 2}, trans
     return μ_pred, σ2_pred
 end
 
-
 """
 svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, Nothing}) where {FT}
 
@@ -524,11 +533,11 @@ function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, No
             println(k)
             # Apply truncated SVD
             n = size(obs_noise_cov)[1]
-            decomposition.S[k+1:n] = zeros(n-k)
             sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S))
-            transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
-            transformed_data = transformed_data[1:k, :];
-            #transformed_data = convert(Matrix{Float64},transformed_data');
+	    transformed_data = sqrt_singular_values_inv[1:k,1:k] * decomposition.Vt[1:k,:] * data
+            transformed_data = transformed_data;
+            decomposition = decomp_struct(decomposition.V[:,1:k], decomposition.Vt[1:k,:], 
+                                   decomposition.S[1:k], n)
 	else
             decomposition = svd(obs_noise_cov)
             sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
@@ -559,11 +568,11 @@ function svd_transform(data::Vector{FT},
             println(k)
             # Apply truncated SVD
             n = size(obs_noise_cov)[1]
-            decomposition.S[k+1:n] = zeros(n-k)
             sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S))
-            transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
-            transformed_data = transformed_data[1:k];
-	    #transformed_data = permutedims(transformed_data, [2, 1]);
+	    transformed_data = sqrt_singular_values_inv[1:k,1:k] * decomposition.Vt[1:k,:] * data
+            transformed_data = transformed_data;
+            decomposition = decomp_struct(decomposition.V[:,1:k], decomposition.Vt[1:k,:], 
+                                   decomposition.S[1:k], n)
 	else
             decomposition = svd(obs_noise_cov)
             sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
@@ -592,9 +601,8 @@ covariance at each point, as a vector of length N_predicted_points, where
 each element is a matrix of size output_dim × output_dim
 """
 function svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::Array{FT, 2}, 
-                                        decomposition::SVD;
+                                        decomposition::Union{SVD, decomp_struct};
 					truncate_svd::FT=1.0) where {FT}
-    if truncate_svd==1.0
         output_dim, N_predicted_points = size(σ2)
         # We created meanvGP = D_inv * Vt * mean_v so meanv = V * D * meanvGP
         sqrt_singular_values= Diagonal(sqrt.(decomposition.S))
@@ -607,9 +615,6 @@ function svd_reverse_transform_mean_cov(μ::Array{FT, 2}, σ2::Array{FT, 2},
             σ2_j = decomposition.V * sqrt_singular_values * Diagonal(σ2[:,j]) * sqrt_singular_values * decomposition.Vt
             transformed_σ2[j] = σ2_j
         end
-    else
-        println("To do")
-    end
 
     return transformed_μ, transformed_σ2
 end

From d059f5f2abb2ba8e5bd22db93845596c305d726c Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Wed, 14 Apr 2021 18:45:45 -0700
Subject: [PATCH 09/10] added tests and lorenz doc

---
 docs/make.jl                             |  4 ++-
 docs/src/examples/lorenz_example.md      | 26 +++++++++++++++++++
 examples/Lorenz/Lorenz_example.jl        | 10 +++++---
 src/GaussianProcessEmulator.jl           |  5 ++--
 src/MarkovChainMonteCarlo.jl             |  4 +--
 test/GaussianProcessEmulator/runtests.jl | 32 +++++++++++++++++++++++-
 test/MarkovChainMonteCarlo/runtests.jl   | 28 ++++++++++++++++++---
 7 files changed, 94 insertions(+), 15 deletions(-)
 create mode 100644 docs/src/examples/lorenz_example.md

diff --git a/docs/make.jl b/docs/make.jl
index b5b8885f6..e5806af48 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -12,11 +12,13 @@ api = ["CalibrateEmulateSample" => [
 				 "Utilities" => "API/Utilities.md"
          ],
 ]
-    
+
+examples = ["Lorenz example" => "examples/lorenz_example.md"]
 
 pages = [
     "Home" => "index.md",
     "Installation instructions" => "installation_instructions.md",
+    "Examples" => examples,
     "API" => api,
 ]
 
diff --git a/docs/src/examples/lorenz_example.md b/docs/src/examples/lorenz_example.md
new file mode 100644
index 000000000..bcf8689d7
--- /dev/null
+++ b/docs/src/examples/lorenz_example.md
@@ -0,0 +1,26 @@
+# Lorenz 96 example
+
+We provide the following template for how the tools may be applied.
+
+For small examples typically have 2 files.
+
+- `GModel.jl` Contains the forward map. The inputs should be the so-called free parameters we are interested in learning, and the output should be the measured data
+- The example script which contains the inverse problem setup and solve
+
+## The structure of the example script
+First we create the data and the setting for the model
+1. Set up the forward model.
+2. Construct/load the truth data. Store this data conveniently in the `Observations.Obs` object
+
+Then we set up the inverse problem
+3. Define the prior distributions. Use the `ParameterDistribution` object
+4. Decide on which `process` tool you would like to use (we recommend you begin with `Invesion()`). Then initialize this with the relevant constructor
+5. initialize the `EnsembleKalmanProcess` object
+
+Then we solve the inverse problem, in a loop perform the following for as many iterations as required:
+7. Obtain the current parameter ensemble
+8. Transform them from the unbounded computational space to the physical space
+9. call the forward map on the ensemble of parameters, producing an ensemble of measured data
+10. call the `update_ensemble!` function to generate a new parameter ensemble based on the new data
+
+One can then obtain the solution, dependent on the `process` type.
diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl
index 9404b26da..1e4e759f0 100644
--- a/examples/Lorenz/Lorenz_example.jl
+++ b/examples/Lorenz/Lorenz_example.jl
@@ -348,8 +348,9 @@ step = 0.1 # first guess
 max_iter = 2000
 yt_sample = truth_sample
 mcmc_test = MarkovChainMonteCarlo.MCMC(yt_sample, Γy, priors, step, u0, max_iter, 
-                         mcmc_alg, burnin, norm_factor;
-                         svdflag=svd_flag, standardize=standardize, truncate_svd=truncate_svd)
+                         mcmc_alg, burnin;
+                         svdflag=svd_flag, standardize=standardize, truncate_svd=truncate_svd,
+			 norm_factor=norm_factor)
 new_step = MarkovChainMonteCarlo.find_mcmc_step!(mcmc_test, gpobj, max_iter=max_iter)
 
 # Now begin the actual MCMC
@@ -360,8 +361,9 @@ burnin = 2000
 max_iter = 100000
 
 mcmc = MarkovChainMonteCarlo.MCMC(yt_sample, Γy, priors, new_step, u0, max_iter, 
-                    mcmc_alg, burnin, norm_factor;
-                    svdflag=svd_flag, standardize=standardize, truncate_svd=truncate_svd)
+                    mcmc_alg, burnin;
+                    svdflag=svd_flag, standardize=standardize, truncate_svd=truncate_svd,
+		    norm_factor=norm_factor)
 MarkovChainMonteCarlo.sample_posterior!(mcmc, gpobj, max_iter)
 
 println("Posterior size")
diff --git a/src/GaussianProcessEmulator.jl b/src/GaussianProcessEmulator.jl
index ce9d8545a..b0d1c511e 100644
--- a/src/GaussianProcessEmulator.jl
+++ b/src/GaussianProcessEmulator.jl
@@ -516,11 +516,10 @@ Note: If F::SVD is the factorization object, U, S, V and Vt can be obtained via
 F.U, F.S, F.V and F.Vt, such that A = U * Diagonal(S) * Vt. The singular values 
 in S are sorted in descending order.
 """
-# TO DO: Add truncate_svd as an input flag here
 function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, Nothing}; 
 		       truncate_svd::FT=1.0) where {FT}
     if obs_noise_cov != nothing
-        # MFH, 3/22/21: Truncate the SVD as a form of regularization
+        # Truncate the SVD as a form of regularization
 	if truncate_svd<1.0 # this variable needs to be provided to this function
             # Perform SVD
             decomposition = svd(obs_noise_cov)
@@ -555,7 +554,7 @@ function svd_transform(data::Vector{FT},
 		       obs_noise_cov::Union{Array{FT, 2}, Nothing};
 		       truncate_svd::FT=1.0) where {FT}
      if obs_noise_cov != nothing
-        # MFH, 3/22/21: Truncate the SVD as a form of regularization
+        # Truncate the SVD as a form of regularization
 	if truncate_svd<1.0 # this variable needs to be provided to this function
             # Perform SVD
             decomposition = svd(obs_noise_cov)
diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl
index d8018eede..64c54f7d1 100644
--- a/src/MarkovChainMonteCarlo.jl
+++ b/src/MarkovChainMonteCarlo.jl
@@ -74,10 +74,10 @@ function MCMC(
     param_init::Vector{FT},
     max_iter::IT,
     algtype::String,
-    burnin::IT,
-    norm_factor::Union{Array{FT, 1}, Nothing};
+    burnin::IT;
     svdflag=true,
     standardize=false,
+    norm_factor::Union{Array{FT, 1}, Nothing}=nothing,
     truncate_svd=1.0) where {FT<:AbstractFloat, IT<:Int}
 
     
diff --git a/test/GaussianProcessEmulator/runtests.jl b/test/GaussianProcessEmulator/runtests.jl
index 64ae66bfc..9b2bf4390 100644
--- a/test/GaussianProcessEmulator/runtests.jl
+++ b/test/GaussianProcessEmulator/runtests.jl
@@ -5,6 +5,7 @@ using GaussianProcesses
 using Statistics 
 using Distributions
 using ScikitLearn
+using LinearAlgebra
 @sk_import gaussian_process : GaussianProcessRegressor
 @sk_import gaussian_process.kernels : (RBF, WhiteKernel, ConstantKernel)
 
@@ -70,6 +71,8 @@ using CalibrateEmulateSample.DataStorage
     @test gp1_norm.sqrt_inv_input_cov ≈ [sqrt(1.0 / Statistics.var(x))] atol=1e-4
 
 
+        
+    
     # GaussianProcess 2: GPJL, predict_f
     pred_type = FType()
     gp2 = GaussianProcess(iopairs, gppackage; GPkernel=GPkernel, obs_noise_cov=nothing, 
@@ -97,7 +100,7 @@ using CalibrateEmulateSample.DataStorage
     μ3, σ3² = GaussianProcessEmulator.predict(gp3, new_inputs)
     @test vec(μ3) ≈ [0.0, 1.0, 0.0, -1.0, 0.0] atol=0.3
     @test vec(σ3²) ≈ [0.016, 0.002, 0.003, 0.004, 0.003] atol=1e-2
-    
+  
 
     # -------------------------------------------------------------------------
     # Test case 2: 2D input, 2D output
@@ -156,5 +159,32 @@ using CalibrateEmulateSample.DataStorage
     @test μ4[:, 4] ≈ [0.0, -2.0] atol=0.25
     @test length(σ4²) == size(new_inputs,2)
     @test size(σ4²[1]) == (d, d)
+
+    # Check if standardization works
+    norm_factor = 10.0
+    norm_factor = fill(norm_factor, size(Y[:,1])) # must be size of output dim
+    println(norm_factor)
+    gp4_standardized = GaussianProcess(iopairs2, gppackage, GPkernel=nothing, obs_noise_cov=Σ, 
+                          normalized=true, noise_learn=true, 
+		          truncate_svd=1.0, standardize=true,
+                          prediction_type=pred_type, norm_factor=norm_factor)
+    cov_est = gp4_standardized.decomposition.V * diagm(gp4_standardized.decomposition.S) *
+            gp4_standardized.decomposition.Vt
+    @test cov_est ≈ Σ ./ (norm_factor .* norm_factor') atol=1e-2
+    @test cov_est ≈ Σ ./ 100.0 atol=1e-2
+
+    # Check if truncation works
+    norm_factor = 10.0
+    norm_factor = fill(norm_factor, size(Y[:,1])) # must be size of output dim
+    println(norm_factor)
+    gp4_trunc = GaussianProcess(iopairs2, gppackage, GPkernel=nothing, obs_noise_cov=Σ, 
+                          normalized=true, noise_learn=true, 
+		          truncate_svd=0.1, standardize=true,
+                          prediction_type=pred_type, norm_factor=norm_factor)
+    μ4_trunc, σ4_trunc = GaussianProcessEmulator.predict(gp4_trunc, new_inputs, transform_to_real=false)
+    μ4_trunc_real, σ4_trunc_real = GaussianProcessEmulator.predict(gp4_trunc, new_inputs, transform_to_real=true)
+    @test size(μ4_trunc,1) == 1
+    @test size(μ4_trunc_real,1) == 2
     
+
 end
diff --git a/test/MarkovChainMonteCarlo/runtests.jl b/test/MarkovChainMonteCarlo/runtests.jl
index a9d66f245..dbd0e8b67 100644
--- a/test/MarkovChainMonteCarlo/runtests.jl
+++ b/test/MarkovChainMonteCarlo/runtests.jl
@@ -60,8 +60,8 @@ using CalibrateEmulateSample.DataStorage
     max_iter = 5000
     norm_factors=nothing
     mcmc_test = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, 
-                        mcmc_alg, burnin, norm_factors; svdflag=true, standardize=false,
-			truncate_svd=1.0)
+                        mcmc_alg, burnin; svdflag=true, standardize=false,
+			truncate_svd=1.0, norm_factor=norm_factors)
     new_step = find_mcmc_step!(mcmc_test, gp)
 
     # reset parameters 
@@ -70,8 +70,8 @@ using CalibrateEmulateSample.DataStorage
 
     # Now begin the actual MCMC
     mcmc = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, 
-                   mcmc_alg, burnin, norm_factors; svdflag=true, standardize=false,
-                   truncate_svd=1.0)
+                   mcmc_alg, burnin; svdflag=true, standardize=false,
+                   truncate_svd=1.0, norm_factor=norm_factors)
     sample_posterior!(mcmc, gp, max_iter)
     posterior_distribution = get_posterior(mcmc)      
     #post_mean = mean(posterior, dims=1)[1]
@@ -89,4 +89,24 @@ using CalibrateEmulateSample.DataStorage
                                    max_iter, "gibbs", burnin)
     @test isapprox(posterior_mean[1] - π/2, 0.0; atol=4e-1)
 
+    # Standardization and truncation
+    norm_factor = 10.0
+    norm_factor = fill(norm_factor, size(y[:,1])) # must be size of output dim
+    gp = GaussianProcess(iopairs, gppackage; GPkernel=GPkernel, obs_noise_cov=σ2_y, 
+                  normalized=false, noise_learn=true, standardize=true, truncate_svd=0.9, 
+                  prediction_type=pred_type, norm_factor=norm_factor) 
+    mcmc_test = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, 
+                        mcmc_alg, burnin;
+			svdflag=true, standardize=true, norm_factor=norm_factor, truncate_svd=0.9)
+
+    # Now begin the actual MCMC
+    mcmc = MCMC(obs_sample, σ2_y, prior, step, param_init, max_iter, 
+                   mcmc_alg, burnin; svdflag=true, standardize=false,
+                   truncate_svd=1.0, norm_factor=norm_factor)
+    sample_posterior!(mcmc, gp, max_iter)
+    posterior_distribution = get_posterior(mcmc)      
+    #post_mean = mean(posterior, dims=1)[1]
+    posterior_mean2 = get_mean(posterior_distribution)
+    @test posterior_mean2 ≈ posterior_mean atol=0.1
+
 end

From abc2decf07d73535eefccd060180d032990d6f21 Mon Sep 17 00:00:00 2001
From: mhowlan3 <mike.howland13@gmail.com>
Date: Thu, 15 Apr 2021 13:41:51 -0700
Subject: [PATCH 10/10] converted to decomp_struct

---
 src/GaussianProcessEmulator.jl | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/GaussianProcessEmulator.jl b/src/GaussianProcessEmulator.jl
index b0d1c511e..ab567b1b9 100644
--- a/src/GaussianProcessEmulator.jl
+++ b/src/GaussianProcessEmulator.jl
@@ -22,6 +22,7 @@ export predict
 export GPJL, SKLJL
 export YType, FType
 export svd_transform, svd_reverse_transform_mean_cov
+export decomp_struct
 
 """
     GaussianProcessesPackage
@@ -47,14 +48,18 @@ struct YType <: PredictionType end
 struct FType <: PredictionType end
 
 
-# SVD decompsotion structure
+# SVD decomposition structure
 struct decomp_struct{FT<:AbstractFloat, IT<:Int}
     V::Array{FT,2}
     Vt::Array{FT,2}
     S::Array{FT}
     N::IT
 end
-
+# SVD decomposition constructor
+function decomp_struct(svd::SVD)
+	# svd.V is of type adjoint, transformed to Array with [:,:]
+	return decomp_struct(svd.V[:,:], svd.Vt, svd.S, size(svd.S)[1])
+end
 
 """
     GaussianProcess{FT<:AbstractFloat}
@@ -76,7 +81,7 @@ struct GaussianProcess{FT<:AbstractFloat, GPM}
     "the Gaussian Process (GP) Regression model(s) that are fitted to the given input-data pairs"
     models::Vector
     "the singular value decomposition of obs_noise_cov, such that obs_noise_cov = decomposition.U * Diagonal(decomposition.S) * decomposition.Vt."
-    decomposition::Union{SVD, decomp_struct, Nothing}
+    decomposition::Union{decomp_struct, Nothing}
     "whether to fit GP models on normalized inputs ((inputs - input_mean) * sqrt_inv_input_cov)"
     normalized::Bool
     "prediction type (`y` to predict the data, `f` to predict the latent function)"
@@ -118,14 +123,10 @@ function GaussianProcess(
     models = Any[]
    
 
-    # TO DO: Standardize the data here
-    # Can use the full time median or some user define function?
+    # Can use the full time median or some user define function
     if standardize
-        #norm_factors = get_standarizing_factors() 
-	println(size(get_outputs(input_output_pairs)))
 	output_values = get_outputs(input_output_pairs) ./ norm_factor
 	obs_noise_cov = obs_noise_cov ./ (norm_factor .* norm_factor')
-	println(size(output_values))
     else
 	output_values = get_outputs(input_output_pairs)
     end
@@ -146,8 +147,6 @@ function GaussianProcess(
     # (if obs_noise_cov==nothing, transformed_data is equal to data)
     transformed_data, decomposition = svd_transform(output_values, 
 						    obs_noise_cov, truncate_svd=truncate_svd)
-    println(size(transformed_data))
-    #println(transformed_data)
 
     # Overwrite the input-output pairs because of the size discrepency
     if truncate_svd<1.0
@@ -268,8 +267,7 @@ function GaussianProcess(
     # Number of models (We are fitting one model per output dimension)
     N_models = output_dim
 
-    # TO DO: Standardize the data here
-    # Can use the full time median or some user define function?
+    # Can use the full time median or some user define function
     if standardize
         #norm_factors = get_standarizing_factors() 
 	output_values = get_outputs(input_output_pairs) ./ norm_factor
@@ -541,6 +539,7 @@ function svd_transform(data::Array{FT, 2}, obs_noise_cov::Union{Array{FT, 2}, No
             decomposition = svd(obs_noise_cov)
             sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
             transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
+	    decomposition = decomp_struct(svd(obs_noise_cov))
         end
     else
         decomposition = nothing
@@ -576,6 +575,7 @@ function svd_transform(data::Vector{FT},
             decomposition = svd(obs_noise_cov)
             sqrt_singular_values_inv = Diagonal(1.0 ./ sqrt.(decomposition.S)) 
             transformed_data = sqrt_singular_values_inv * decomposition.Vt * data
+	    decomposition = decomp_struct(svd(obs_noise_cov))
         end
     else
         decomposition = nothing