truncated SVD L96 example

.gitignore

@@ -14,3 +14,5 @@ docs/site/
 *.jl.mem
 *.vscode*
 *.DS_Store
+*.jld2
+output*

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

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 = "ω"

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"

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 

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)

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