format nb

notebooks/Linear system with EnKF.ipynb

@@ -2,28 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "\"/media/mat/HDD/EnKF/Project.toml\""
-      ]
-     },
-     "execution_count": 1,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "using Pkg\n",
-    "Pkg.activate(\"/media/mat/HDD/EnKF/\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 228,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -34,7 +13,7 @@
     "using LinearAlgebra\n",
     "using Plots\n",
     "using ProgressMeter\n",
-    "using DifferentialEquations"
+    "using OrdinaryDiffEq"
    ]
   },
   {

notebooks/Lorenz attractor with EnKF RTPS.ipynb

@@ -33,7 +33,7 @@
     "using DocStringExtensions\n",
     "using LinearAlgebra\n",
     "using ProgressMeter\n",
-    "using DifferentialEquations"
+    "using OrdinaryDiffEq"
    ]
   },
   {

notebooks/Lorenz attractor with EnKF.ipynb

@@ -38,7 +38,7 @@
       "└ @ Base.Docs docs/Docs.jl:223\n",
       "┌ Warning: Replacing docs for `EnKF.A :: Union{Tuple{EnsembleState{N,TS}}, Tuple{TS}, Tuple{NS}, Tuple{N}} where TS where NS where N` in module `EnKF`\n",
       "└ @ Base.Docs docs/Docs.jl:223\n",
-      "┌ Info: Recompiling stale cache file /home/mat/.julia/compiled/v1.1/DifferentialEquations/UQdwS.ji for DifferentialEquations [0c46a032-eb83-5123-abaf-570d42b7fbaa]\n",
+      "┌ Info: Recompiling stale cache file /home/mat/.julia/compiled/v1.1/OrdinaryDiffEq/UQdwS.ji for OrdinaryDiffEq [0c46a032-eb83-5123-abaf-570d42b7fbaa]\n",
       "└ @ Base loading.jl:1184\n"
      ]
     }
@@ -50,7 +50,7 @@
     "using DocStringExtensions\n",
     "using LinearAlgebra\n",
     "using ProgressMeter\n",
-    "using DifferentialEquations"
+    "using OrdinaryDiffEq"
    ]
   },
   {

notebooks/Lorenz95 with EnKF.ipynb

@@ -38,7 +38,7 @@
       "└ @ Base.Docs docs/Docs.jl:223\n",
       "┌ Warning: Replacing docs for `EnKF.A :: Union{Tuple{EnsembleState{N,TS}}, Tuple{TS}, Tuple{NS}, Tuple{N}} where TS where NS where N` in module `EnKF`\n",
       "└ @ Base.Docs docs/Docs.jl:223\n",
-      "┌ Info: Recompiling stale cache file /home/mat/.julia/compiled/v1.1/DifferentialEquations/UQdwS.ji for DifferentialEquations [0c46a032-eb83-5123-abaf-570d42b7fbaa]\n",
+      "┌ Info: Recompiling stale cache file /home/mat/.julia/compiled/v1.1/OrdinaryDiffEq/UQdwS.ji for OrdinaryDiffEq [0c46a032-eb83-5123-abaf-570d42b7fbaa]\n",
       "└ @ Base loading.jl:1184\n"
      ]
     }
@@ -50,7 +50,7 @@
     "using DocStringExtensions\n",
     "using LinearAlgebra\n",
     "using ProgressMeter\n",
-    "using DifferentialEquations"
+    "using OrdinaryDiffEq"
    ]
   },
   {

notebooks/Propagator using DifferentialEquations.ipynb

@@ -42,7 +42,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "using DifferentialEquations"
+    "using OrdinaryDiffEq"
    ]
   },
   {

notebooks/Test ETKF.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {

notebooks/VortexEnKF.jl

@@ -9,11 +9,9 @@ using PotentialFlow
 """
 Create a specific structure for EnKF-based Vortex model,
 Darakananda, et al. Data-assimilated low-order vortex modeling of separated flows, APS (2018)
-
 """
 
 
-
 """
     VortexENKFprop
 
@@ -114,8 +112,6 @@ VortexENKFfilter = FilteringFunction()
 
 
 function VortexENKF(N, NZ, ri, m, z, ϵ)
-
 #     return ENKF{N, NZ}(VortexENKFprop, ri, VortexENKFfilter, m, z, ϵ, true, true, true)
-
     return ENKF{N, NZ}(VortexENKFprop, ri, VortexENKFfilter, m, z, ϵ, true, true, true, true);
 end

src/EnKF.jl

@@ -11,6 +11,7 @@ export PropagationFunction,
   MeasurementFunction,
   FilteringFunction,
   RealMeasurementFunction
+export ENKF
 
 struct PropagationFunction end
 struct MeasurementFunction end
@@ -24,7 +25,185 @@ include("inflation.jl")
 
 # include("update.jl")
 # include("system.jl")
-include("stochasticEnKF.jl")
+# include("stochasticEnKF.jl")
 include("ETKF.jl")
 
+
+""""
+Define system ENKF which performs the data assimilation
+using the stochastic ensemble Kalman filter
+
+Fields:
+ - 'f' : propagation function
+ - 'A' : inflation
+ - 'G' : filtering function acting on the state
+ 
+ - 'm' : measurement function based on state
+ - 'z' : real measurement function
+ - 'ϵ' : measurement noise distribution
+
+ - 'bounds' : bounds on certain states
+
+ - 'isinflated' : Bool = true if state is inflated,
+     = false otherwise
+
+ - 'isfiltered' : Bool = true if state has to be filtered,
+     = false otherwise
+
+ - 'isaugmented' : Bool = true if measurement function is nonlinear,
+     = false otherwise
+"""
+mutable struct ENKF{N,NZ}
+  # "Ensemble of states"
+  # ENS::EnsembleState{N, NS, TS}
+  "Propagation function"
+  f::PropagationFunction
+
+  "Covariance Inflation"
+  A::Union{InflationType,RecipeInflation}
+
+  "Filter function"
+  G::FilteringFunction
+
+  "Measurement function based on state"
+  m::MeasurementFunction
+
+  "Real measurement function"
+  z::RealMeasurementFunction
+
+  "Measurement noise distribution"
+  ϵ::AdditiveInflation{NZ}
+
+  "Boolean: is data assimilation used"
+  isenkf::Bool
+
+  "Boolean: is state vector inflated"
+  isinflated::Bool
+
+  "Boolean: is state vector filtered"
+  isfiltered::Bool
+
+  "Boolean: is state vector augmented"
+  isaugmented::Bool
+end
+# "Bounds on certain state"
+# bounds
+
+
+"""
+    Define action of ENKF on EnsembleState
+"""
+function (enkf::ENKF{N,NZ})(t::Float64,
+  Δt::Float64,
+  ens::EnsembleState{N,TS}) where {N,NZ,TS}
+
+  "Propagate each ensemble member"
+  enkf.f(t, ens)
+
+  "Is data assimilation used"
+  if enkf.isenkf == true
+    # println("good prop")
+    "Covariance inflation if 'isinflated==true' "
+    if enkf.isinflated == true
+      enkf.A(ens)
+    end
+
+    "State filtering if 'isfiltered==true' "
+    if enkf.isfiltered == true
+      enkf.G(ens)
+    end
+    # println("good filtering")
+
+    "Compute mean and deviation"
+    ensfluc = EnsembleState(N, ens.S[1])
+    deviation(ensfluc, ens)
+
+    A′ = hcat(ensfluc)
+
+    " Additional computing for RTPS inflation"
+    if typeof(enkf.A) <: Union{RTPSInflation,RTPSAdditiveInflation,RTPSRecipeInflation}
+      # correct scaling by 1/N-1 instead of 1/N for small ensembles
+      σᵇ = std(ensfluc.S, corrected=false)
+    end
+
+    # println("good deviation")
+    "Compute measurement"
+    mens = EnsembleState(N, zeros(NZ))
+
+    for (i, s) in enumerate(ens.S)
+      mens.S[i] = enkf.m(t, deepcopy(s))
+    end
+    Â = hcat(deepcopy(mens))
+
+    "Define measurement matrix H for linear observations, can be time-varying"
+    if enkf.isaugmented == false
+      H = deepcopy(enkf.m(t))
+    end
+
+    "Compute deviation from measurement of the mean"
+    Ŝ = mean(deepcopy(ens))
+
+    Â′ = Â .- enkf.m(t, Ŝ)
+
+    # println("good deviation mean")
+    "Get actual measurement"
+    zens = EnsembleState(N, zeros(NZ))
+    enkf.z(t + Δt, zens)
+    # @show zens
+    # println("good actual measurement")
+
+    "Perturb actual measurement"
+    enkf.ϵ(zens)
+    D = hcat(zens)
+
+    "Analysis step with representers, Evensen, Leeuwen et al. 1998"
+    "Construct representers"
+    if enkf.isaugmented == true
+      b = ((Â′ * Â′') + (N - 1) * cov(enkf.ϵ) * I) \ (D - Â)
+      Bᵀb = (A′ * Â′') * b
+    else
+      b = (H * (A′ * A′') * H' + (N - 1) * cov(enkf.ϵ) * I) \ (D - Â)
+      Bᵀb = (A′ * A′') * H' * b
+    end
+
+    "Analysis step"
+    ens += cut(Bᵀb)
+
+    " Additional computing for RTPS inflation"
+    if typeof(enkf.A) <: Union{RTPSInflation,RTPSAdditiveInflation,RTPSRecipeInflation}
+      ensfluc = EnsembleState(N, ens.S[1])
+      deviation(ensfluc, ens)
+      σᵃ = std(ensfluc.S, corrected=false)
+      enkf.A(ens, σᵇ, σᵃ)
+    end
+
+    "State filtering if 'isfiltered==true' "
+    if enkf.isfiltered == true
+      enkf.G(ens)
+    end
+
+    " Compute a posteriori covariance"
+    deviation(ensfluc, ens)
+
+    A′ = hcat(ensfluc)
+    return t + Δt, ens, A′ * A′'
+  else
+    return t + Δt, ens
+  end
+end
+
+# Create constructor for ENKF
+function ENKF(N, NZ, f, A, G, m, z, ϵ;
+  isenkf::Bool=true, isinflated::Bool=false, isfiltered::Bool=false, isaugmented::Bool=false)
+  return ENKF{N,NZ}(f, A, G, m, z, ϵ, isenkf, isinflated, isfiltered, isaugmented)
+end
+
+
+# size(enkf::ENKF{N, TS, NZ, TZ}) where {N, TS, NZ, TZ} = (N, size, NZ)
+
+# function Base.show(io::IO, sys::ENKF{N, TS, NZ, TZ}) where {N, TS, NZ, TZ}
+#     NS = size()
+#     print(io, "Ensemble Kalman filter with $N members of state of size $ and measurement vector of length $NZ")
+# end
+
 end # module