diff --git a/Project.toml b/Project.toml
index 4c2c95b47..066b0eaf5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,9 +9,14 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 GaussianProcesses = "891a1506-143c-57d2-908e-e1f8e92e6de9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+
+[targets]
+test = ["Test", "Plots"]
diff --git a/src/EKP.jl b/src/EKP.jl
index e6fd2ca59..72042bd5b 100644
--- a/src/EKP.jl
+++ b/src/EKP.jl
@@ -72,14 +72,6 @@ function EKObj(parameters::Array{FT, 2},
                process::P;
                Δt=FT(1)) where {FT<:AbstractFloat, P<:Process}
 
-    # Throw an error if an attempt is made to instantiate a `Sampler` EKObj
-    # EK Sampler implementation is not finalized yet, so its use is prohibited
-    # TODO: Finalize EKS implementation (can be done as soon as we know the
-    #       correct EKS update equation, which apparently is different from
-    #       Eq. (2.8) in Cleary et al. (2019))
-    # err_msg = "Ensemble Kalman Sampler is not fully implemented yet. Use Ensemble Kalman Inversion instead."
-    # typeof(process) != Sampler{FT} || error(err_msg)
-
     # ensemble size
     N_ens = size(parameters)[1]
     IT = typeof(N_ens)
@@ -252,11 +244,6 @@ function update_ensemble!(ek::EKObj{FT, IT, Sampler{FT}}, g) where {FT, IT}
 
     noise = MvNormal(u_cov)
 
-
-    ###########################################################################
-    ###############    TODO: Implement correct equation here   ################
-    ###########################################################################
-
     implicit = (1 * Matrix(I, size(u)[2], size(u)[2]) + Δt * (ek.process.prior_cov' \ u_cov')') \
                   (u'
                     .- Δt * ( u' .- u_mean) * D
@@ -265,9 +252,6 @@ function update_ensemble!(ek::EKObj{FT, IT, Sampler{FT}}, g) where {FT, IT}
 
     u = implicit' + sqrt(2*Δt) * rand(noise, N_ens)'
 
-    ###########################################################################
-    ###########################################################################
-
     # store new parameters (and observations)
     push!(ek.u, u) # N_ens x N_params
     push!(ek.g, g) # N_ens x N_data