dapper.mods.Lorenz96.sakov2008

  1"""Settings as in `bib.sakov2008deterministic`.
  2
  3This HMM is used (with small variations) in many DA papers, for example
  4
  5`bib.bocquet2011ensemble`, `bib.sakov2012iterative`,
  6`bib.bocquet2015expanding`, `bib.bocquet2013joint`.
  7"""
  8
  9import numpy as np
 10
 11import dapper.mods as modelling
 12from dapper.mods.Lorenz96 import LPs, Tplot, dstep_dx, step, x0
 13from dapper.tools.localization import nd_Id_localization
 14
 15# Sakov uses K=300000, BurnIn=1000*0.05
 16tseq = modelling.Chronology(0.05, dko=1, Ko=1000, Tplot=Tplot, BurnIn=2 * Tplot)
 17
 18Nx = 40
 19x0 = x0(Nx)
 20
 21Dyn = {
 22    "M": Nx,
 23    "model": step,
 24    "linear": dstep_dx,
 25    "noise": 0,
 26}
 27
 28X0 = modelling.GaussRV(mu=x0, C=0.001)
 29
 30jj = np.arange(Nx)  # obs_inds
 31Obs = modelling.partial_Id_Obs(Nx, jj)
 32Obs["noise"] = 1
 33Obs["localizer"] = nd_Id_localization((Nx,), (2,))
 34
 35HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0)
 36
 37HMM.liveplotters = LPs(jj)
 38
 39
 40####################
 41# Suggested tuning
 42####################
 43
 44# Reproduce Table1 of sakov2008deterministic        # Expected rmse.a:
 45# --------------------------------------------------------------------------------
 46# xps += EnKF('PertObs'        ,N=40, infl=1.06)               # 0.22
 47# xps += EnKF('DEnKF'          ,N=40, infl=1.01)               # 0.18
 48# xps += EnKF('PertObs'        ,N=28, infl=1.08)               # 0.24
 49# xps += EnKF('Sqrt'           ,N=24, infl=1.013,rot=True)     # 0.18
 50#
 51# Other analysis schemes:
 52# xps += EnKF('Serial'         ,N=28, infl=1.02,rot=True)      # 0.18
 53# xps += EnKF('Serial ESOPS'   ,N=28, infl=1.02)               # 0.18
 54# xps += EnKF('Serial Stoch'   ,N=28, infl=1.08)               # 0.24
 55#
 56# EnKF-N
 57# xps += EnKF_N(N=24,rot=True)                                 # 0.21
 58# xps += EnKF_N(N=24,rot=True,xN=2.0)                          # 0.18
 59#
 60# Baseline methods
 61# xps += Climatology()                                         # 3.6
 62# xps += OptInterp()                                           # 0.95
 63# xps += Var3D(xB=0.02)                                        # 0.41
 64# xps += ExtKF(infl=10)                                        # 0.24
 65
 66# Reproduce LETKF scores from bocquet2011ensemble fig 6:
 67# --------------------------------------------------------------------------------
 68# xps += LETKF(N=6,rot=True,infl=1.05,loc_rad=4,taper='Step')  #
 69# Other localized:
 70# xps += LETKF(         N=7,rot=True,infl=1.04,loc_rad=4)      # 0.22
 71# xps += SL_EAKF(       N=7,rot=True,infl=1.07,loc_rad=6)      # 0.23
 72
 73# Reproduce Table 3 (IEnKF) from sakov2012iterative
 74# --------------------------------------------------------------------------------
 75# HMM.tseq.dko = 12
 76# xps += iEnKS('Sqrt' ,N=25,Lag=1,nIter=10,infl=1.2,rot=1)     # 0.46
 77# xps += iEnKS('Sqrt' ,N=25,Lag=1,nIter=10,xN=2.0  ,rot=1)     # 0.46
 78
 79# Reproduce Fig 3 of Bocquet'2015 "expanding"
 80# --------------------------------------------------------------------------------
 81# xps += EnKF('Sqrt',N=20,rot=True,infl=1.04)                  # 0.20
 82# # use infl=1.10 with dko=3
 83# # use infl=1.40 with dko=5
 84# xps += EnKF_N(N=20)                                          # 0.24
 85# xps += EnKF_N(N=20,xN=2)                                     # 0.19
 86# # Also try quasi-linear regime:
 87# tseq = Chronology(0.01,dko=1,...)
 88
 89# Reproduce Bocquet/Sakov'2013 "Joint...", Fig 4, i.e. dto=0.2:
 90# xps += iEnKS('Sqrt', N=20, Lag=4, xN=2) # 0.31
 91# xps += Var4D(Lag=1,xB=0.2)              # 0.46
 92# xps += Var4D(Lag=2,xB=0.1)              # 0.39
 93# xps += Var4D(Lag=4,xB=0.02)             # 0.37
 94# Cannot reproduce Fig4's reported 4D-Var scores for L>4. Example:
 95# xps += Var4D(Lag=6,xB=0.015)            # 0.385 Boc13 reports 0.33
 96
 97# Tests with the Particle filter, with N=3000, Ko=10'000.
 98# da_method  NER  reg  |  rmse.a   rmv.a
 99# --------- ----  ---  -  ------  ------
100# PartFilt  0.05  1.2  |  0.35    0.40
101# PartFilt  0.05  1.6  |  0.41    0.45
102# PartFilt  0.5   0.7  |  0.26    0.29
103# PartFilt  0.5   0.9  |  0.30    0.34
104# PartFilt  0.5   1.2  |  0.36    0.40
105# Using NER=0.9 yielded rather poor results.