Merge #266

266: Small bug fixes in tests r=odunbar a=odunbar

<!--- THESE LINES ARE COMMENTED -->
## Purpose 
<!--- One sentence to describe the purpose of this PR, refer to any linked issues:
#14 -- this will link to issue 14
Closes #2 -- this will automatically close issue 2 on PR merge
-->
- Closes #265 
- Addresses #261 (1/2)
Leads to a more robust test pipeline and constructor `constrained_gaussian`

## Content
<!---  specific tasks that are currently complete 
- Solution implemented
-->
- fixes bugs within #265 ,
- resolves an eks object bug where it would mix prior and final ensemble, and now puts plots in the EKP test subdirectory
- adds warning for when users prescribe values in `constrained_gaussian` which will be slow to solve. In these cases a large amount of mass will be concentrated on the boundary.
- tightens error of termination to `mu + sig < upper_bound`  (before was `mu + 0.7*sig > 1` etc.) likewise for lower bound

<!---
Review checklist

I have:
- followed the codebase contribution guide: https://clima.github.io/ClimateMachine.jl/latest/Contributing/
- followed the style guide: https://clima.github.io/ClimateMachine.jl/latest/DevDocs/CodeStyle/
- followed the documentation policy: https://github.com/CliMA/policies/wiki/Documentation-Policy
- checked that this PR does not duplicate an open PR.

In the Content, I have included 
- relevant unit tests, and integration tests, 
- appropriate docstrings on all functions, structs, and modules, and included relevant documentation.

-->

----
- [x] I have read and checked the items on the review checklist.


Co-authored-by: odunbar <[email protected]>

src/ParameterDistributions.jl

@@ -979,12 +979,23 @@ function constrained_gaussian(
             μ_u, σ_u = _inverse_lognormal_mean_std(μ_c - lower_bound, σ_c)
         else
             # finite interval case; need to solve numerically
-            if (μ_c - 0.7 * σ_c <= lower_bound)
+            if (μ_c - σ_c <= lower_bound)
                 throw(DomainError("`$(name)`: Target std $(σ_c) puts μ - σ too close to lower bound $(lower_bound)"))
             end
-            if (μ_c + 0.7 * σ_c >= upper_bound)
+            if (μ_c + σ_c >= upper_bound)
                 throw(DomainError("`$(name)`: Target std $(σ_c) puts μ + σ too close to upper bound $(upper_bound)"))
             end
+            # 1.2 seems a reasonable tolerance here for solver to converge quickly
+            if (μ_c - 1.2 * σ_c <= lower_bound)
+                @warn(
+                    "`$(name)`: Target std $(σ_c) puts μ - σ very close to lower bound $(lower_bound), \n The solver may need more iterations to converge, consider decreasing σ"
+                )
+            end
+            if (μ_c + 1.2 * σ_c >= upper_bound)
+                @warn(
+                    "`$(name)`: Target std $(σ_c) puts μ + σ very close to upper bound $(upper_bound), \n The solver may need more iterations to converge, consider decreasing σ"
+                )
+            end
             μ_u, σ_u = _constrained_gaussian(
                 name,
                 μ_c,

test/EnsembleKalmanProcess/runtests.jl

@@ -9,6 +9,8 @@ using EnsembleKalmanProcesses.Localizers
 import EnsembleKalmanProcesses: construct_mean, construct_cov, construct_sigma_ensemble
 const EKP = EnsembleKalmanProcesses
 
+
+
 # Read inverse problem definitions
 include("inverse_problem.jl")
 
@@ -122,10 +124,11 @@ end
     # Plot evolution of the EKS particles
     if TEST_PLOT_OUTPUT
         gr()
+        eksobj = eksobjs[end]
         ϕ_prior = transform_unconstrained_to_constrained(prior, get_u_prior(eksobj))
         ϕ_final = get_ϕ_final(prior, eksobj)
-        p = plot(ϕ_prior[1, :], ϕ_final[2, :], seriestype = :scatter, label = "Initial ensemble")
-        plot!(ϕ_prior[1, :], ϕ_final[2, :], seriestype = :scatter, label = "Final ensemble")
+        p = plot(ϕ_prior[1, :], ϕ_prior[2, :], seriestype = :scatter, label = "Initial ensemble")
+        plot!(ϕ_final[1, :], ϕ_final[2, :], seriestype = :scatter, label = "Final ensemble")
         plot!(
             [ϕ_star[1]],
             xaxis = "cons_p",
@@ -136,7 +139,7 @@ end
             label = :none,
         )
         plot!([ϕ_star[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = :none)
-        savefig(p, "EKS_test.png")
+        savefig(p, joinpath(@__DIR__, "EKS_test.png"))
     end
 end
 
@@ -290,7 +293,7 @@ end
             label = :none,
         )
         plot!([ϕ_star[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = :none)
-        savefig(p, "EKI_test.png")
+        savefig(p, joinpath(@__DIR__, "EKI_test.png"))
     end
 end
 
@@ -393,18 +396,18 @@ end
     if TEST_PLOT_OUTPUT
         gr()
         θ_mean_arr = hcat(ukiobj.process.u_mean...)
-        N_θ = length(ukiobj.process.u_mean[1])
-        θθ_std_arr = zeros(Float64, (N_θ, N_iter + 1))
-        for i in 1:(N_iter + 1)
+        N_θ, N_ens = size(θ_mean_arr)
+        θθ_std_arr = zeros(Float64, (N_θ, N_ens))
+        for i in 1:(N_ens)
             for j in 1:N_θ
                 θθ_std_arr[j, i] = sqrt(ukiobj.process.uu_cov[i][j, j])
             end
         end
 
         u_star = transform_constrained_to_unconstrained(prior, ϕ_star)
-        ites = Array(LinRange(1, N_iter + 1, N_iter + 1))
+        ites = Array(LinRange(1, N_ens, N_ens))
         p = plot(ites, grid = false, θ_mean_arr[1, :], yerror = 3.0 * θθ_std_arr[1, :], label = "cons_p")
-        plot!(ites, fill(u_star[1], N_iter + 1), linestyle = :dash, linecolor = :grey, label = :none)
+        plot!(ites, fill(u_star[1], N_ens), linestyle = :dash, linecolor = :grey, label = :none)
         plot!(
             ites,
             grid = false,
@@ -413,8 +416,8 @@ end
             label = "uncons_p",
             xaxis = "Iterations",
         )
-        plot!(ites, fill(u_star[2], N_iter + 1), linestyle = :dash, linecolor = :grey, label = :none)
-        savefig(p, "UKI_test.png")
+        plot!(ites, fill(u_star[2], N_ens), linestyle = :dash, linecolor = :grey, label = :none)
+        savefig(p, joinpath(@__DIR__, "UKI_test.png"))
     end
 end
 

test/Inflation/runtests.jl

@@ -82,11 +82,11 @@ initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ens)
             EKP.additive_inflation!(eki_add_inflation_prior; use_prior_cov = true)
 
             # ensure multiplicative inflation approximately preserves ensemble mean
-            @test get_u_mean_final(ekiobj) ≈ get_u_mean_final(eki_mult_inflation) atol = 0.1
+            @test get_u_mean_final(ekiobj) ≈ get_u_mean_final(eki_mult_inflation) atol = 0.2
             # ensure additive inflation approximately preserves ensemble mean
-            @test get_u_mean_final(ekiobj) ≈ get_u_mean_final(eki_add_inflation) atol = 0.1
+            @test get_u_mean_final(ekiobj) ≈ get_u_mean_final(eki_add_inflation) atol = 0.2
             # ensure additive inflation (scaling prior cov) approximately preserves ensemble mean
-            @test get_u_mean_final(ekiobj) ≈ get_u_mean_final(eki_add_inflation_prior) atol = 0.1
+            @test get_u_mean_final(ekiobj) ≈ get_u_mean_final(eki_add_inflation_prior) atol = 0.2
 
             # ensure inflation expands ensemble variance as expected
             expected_var_gain = 1 / (1 - Δt)

test/ParameterDistributions/runtests.jl

@@ -617,6 +617,11 @@ using EnsembleKalmanProcesses.ParameterDistributions
             @test_throws DomainError constrained_gaussian("test", 0.0, 10.0, -1.0, 1000.0)
             @test_throws DomainError constrained_gaussian("test", 0.0, 10.0, -1000.0, 1.0)
             @test_throws DomainError constrained_gaussian("test", 0.0, 10.0, -1.0, 1.0)
+            # σ near boundary throws warning
+            @test_logs (:warn,) constrained_gaussian("test", 0.54, 0.4, 0, 1) # 0.54 + 1.2*0.4 > 1
+            @test_logs (:warn,) constrained_gaussian("test", 0.46, 0.4, 0, 1) # 0.46 - 1.2*0.4 < 1
+
+
         end
         @testset "constrained_gaussian: closed form" begin
             μ_c = -5.0