add example

examples/ex_duffing.jl

@@ -0,0 +1,66 @@
+using EnKF
+using OrdinaryDiffEq
+using Plots
+gr(framestyle=:box)
+
+include("main_duffing.jl")
+
+u0 = [1.0; -1.0]
+tspan = (0.0, 50.0)
+
+Δt = 1e-2
+T = tspan[1]:Δt:tspan[end]
+
+prob = ODEProblem(duffing, u0, tspan)
+sol = solve(prob, RK4(), adaptive=false, dt=Δt)
+
+integrator = init(prob, RK4(), adaptive=false, dt=Δt, save_everystep=false)
+
+plot(sol, idxs=(1))
+
+begin
+  fprop = PropagationFunction()
+  m = MeasurementFunction()
+  z = RealMeasurementFunction()
+  A = IdentityInflation()
+  ϵ = AdditiveInflation(MvNormal(zeros(1), 1.0 * I))
+
+  N = 50
+  NZ = 1
+
+  x₀ = [0.5, -0.5]
+  ens = initialize(N, MvNormal(x₀, 2.0 * I))
+  ens.S[1]
+  estimation_state = [deepcopy(ens.S)]
+
+  true_state = [deepcopy(x₀)]
+  covs = []
+  g = FilteringFunction()
+
+  enkf = ENKF(N, NZ, fprop, A, g, m, z, ϵ;
+    isinflated=false, isfiltered=false, isaugmented=false)
+end
+
+
+Δt = 1e-2
+Tsub = 0.0:Δt:50.0-Δt
+
+for (n, t) in enumerate(Tsub)
+  global ens
+  t, ens, cov = enkf(t, Δt, ens)
+
+  push!(estimation_state, deepcopy(ens.S))
+  push!(covs, deepcopy(cov))
+end
+
+s = hcat(sol(T).u...)
+ŝ = hcat(mean.(estimation_state)...)
+
+
+
+plt = plot(layout=(2, 1), legend=:bottomright)
+plot!(plt[1], T, s[1, 1:end], linewidth=3, label="truth")
+plot!(plt[1], Tsub, ŝ[1, 1:end-1], linewidth=3, markersize=2, label="EnKF mean", xlabel="t", ylabel="x", linestyle=:dash)
+
+plot!(plt[2], T, s[2, 1:end], linewidth=3, label="truth")
+plot!(plt[2], Tsub, ŝ[2, 1:end-1], linewidth=3, markersize=2, label="EnKF mean", xlabel="t", ylabel="y", linestyle=:dash)

examples/main_duffing.jl

@@ -0,0 +1,38 @@
+using EnKF
+
+
+function duffing(du, u, p, t)
+  γ = 0.1
+  d = 0.1
+  ω = 1.4
+
+  du[1] = u[2]
+  du[2] = u[1] - u[1]^3 - γ * u[2] + d * cos(ω * t)
+end
+
+
+
+function (::PropagationFunction)(t::Float64, ENS::EnsembleState{N,TS}) where {N,TS}
+  for (i, s) in enumerate(ENS.S)
+    set_t!(integrator, deepcopy(t))
+    set_u!(integrator, deepcopy(s))
+    step!(integrator)
+    ENS.S[i] = deepcopy(integrator.u)
+  end
+  return ENS
+end
+
+function (::MeasurementFunction)(t::Float64, s::TS) where TS
+  return [s[2]]
+end
+
+function (::MeasurementFunction)(t::Float64) # error function
+  return reshape([0.0, 1.0], (1, 2))
+end
+
+function (::RealMeasurementFunction)(t::Float64, ENS::EnsembleState{N,TZ}) where {N,TZ}
+  let s = sol(t)
+    fill!(ENS, [deepcopy(s[2])])
+  end
+  return ENS
+end

src/EnKF.jl

@@ -1,6 +1,10 @@
 module EnKF
 
 using Distributions, Statistics, LinearAlgebra
+export MvNormal, I
+export mean, std
+
+
 import Statistics: mean, var, std, cov
 import Base: size, length, show, hcat, +, -, fill!
 import Distributions