working example

README.md

@@ -44,7 +44,22 @@ What makes EKP different?
 
 ## What does it look like to use?
 
-Below we will outline the current user experience for using `EnsembleKalmanProcesses.jl` to solve the classic inverse problem where we learn `y = G(u) + e`, for `e` distributed with `N(0,Γ)`. We assume some prior knowledge of the parameters `u` in the problem (say we have five, one positive and four more that are more uncertain), then we are ready to go! 
+Below we will outline the current user experience for using `EnsembleKalmanProcesses.jl` to solve the classic inverse problem where we learn `y = G(u)`, noisy forward map `G` distributed as `N(0,Γ)`. For this first example
+```julia
+using LinearAlgebra
+G(u) = [
+    1/abs(u[1]),
+    sum(u[2:5]),
+    prod(u[3:4]),
+    u[1]^2-u[2]-u[3],
+    u[4],
+    u[5]^3,
+    ] .+ 0.1*randn(6)
+true_u = [3, 1, 2,-3,-4]
+y = G(true_u)
+Γ = (0.1)^2*I
+```
+We assume some prior knowledge of the parameters `u` in the problem (say we have five, one positive and four more that are more uncertain), then we are ready to go! 
 
 ```julia
 using EnsembleKalmanProcesses
@@ -53,16 +68,13 @@ using EnsembleKalmanProcesses.ParameterDistributions
 prior_u1 = constrained_gaussian("positive_with_mean_2", 2, 1, 0, Inf)
 prior_u2 = constrained_gaussian("four_with_spread_5", 0, 5, -Inf, Inf, repeats=4)
 prior = combine_distributions([prior_u1, prior_u2]) 
-using Plots
-plot(prior) # lets see it with Plots.jl
 
-N_ensemble = 10
+N_ensemble = 50
 initial_ensemble = construct_initial_ensemble(prior, N_ensemble)
 ensemble_kalman_process = EnsembleKalmanProcess(
-    initial_ensemble, y, Γ, Inversion() # use Ensemble Kalman Inversion updates
-)
+    initial_ensemble, y, Γ, Inversion(), verbose=true)
 
-N_iterations = 5
+N_iterations = 10
 for i in 1:N_iterations
     params_i = get_ϕ_final(prior, ensemble_kalman_process)
 
@@ -73,11 +85,20 @@ for i in 1:N_iterations
     update_ensemble!(ensemble_kalman_process, G_matrix)
 end
 
-final_solution = get_ϕ_mean_final(prior, ensemble_kalman_process) 
+final_solution = get_ϕ_mean_final(prior, ensemble_kalman_process)
+
+
+# Let's see what's going on!
+using Plots
+p = plot(prior)
+for (i,sp) in enumerate(p.subplots)
+    vline!(sp, [true_u[i]], lc="black", lw=4)
+    vline!(sp, [final_solution[i]], lc="magenta", lw=4)
+end
+display(p)
 ```
-With suitable parallelization, the final solution is obtained in the time it takes to perform ~5 iterations, `G` and five EKP update steps.
 
-See this example working [here!](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/literated/sinusoid_example/). check out our many example scripts above in `examples/`
+See a similar working example [here!](https://clima.github.io/EnsembleKalmanProcesses.jl/dev/literated/sinusoid_example/). Check out our many example scripts above in `examples/`
 
 # [Quick links!](@id quick-links)