dapper.mods.QG.sakov2008

  1"""Reproduce results from `bib.sakov2008deterministic`."""
  2
  3import numpy as np
  4
  5import dapper.mods as modelling
  6from dapper.mods.QG import LP_setup, model_config, sample_filename, shape
  7from dapper.tools.localization import nd_Id_localization
  8
  9############################
 10# Time series, model, initial condition
 11############################
 12
 13model = model_config("sakov2008", {})
 14Dyn = {
 15    "M": np.prod(shape),
 16    "model": model.step,
 17    "noise": 0,
 18}
 19
 20# Considering that I have 8GB mem on the Mac, and the estimate:
 21# ≈ (8 bytes/float)*(129² float/stat)*(7 stat/k) * K,
 22# it should be possible to run experiments of length (K) < 8000.
 23tseq = modelling.Chronology(dt=model.prms["dtout"], dko=1, T=1500, BurnIn=250)
 24# In my opinion the burn in should be 400.
 25# Sakov also used 10 repetitions.
 26
 27X0 = modelling.RV(M=Dyn["M"], file=sample_filename)
 28
 29
 30############################
 31# Observation settings
 32############################
 33
 34# This will look like satellite tracks when plotted in 2D
 35Ny = 300
 36jj = modelling.linspace_int(Dyn["M"], Ny)
 37
 38# Want: random_offset(t1)==random_offset(t2) if t1==t2.
 39# Solutions: (1) use caching (ensure maxsize=inf) or (2) stream seeding.
 40# Either way, use a local random stream to avoid interfering with global stream
 41# (and e.g. ensure equal outcomes for 1st and 2nd run of the python session).
 42# TODO 5: upgrade to `default_rng()`
 43rstream = np.random.RandomState()
 44max_offset = jj[1] - jj[0]
 45
 46
 47def obs_inds(ko):
 48    def random_offset():
 49        rstream.seed(ko)
 50        u = rstream.rand()
 51        return int(np.floor(max_offset * u))
 52
 53    return jj + random_offset()
 54
 55
 56def obs_now(ko):
 57    jj = obs_inds(ko)
 58
 59    @modelling.ens_compatible
 60    def hmod(E):
 61        return E[jj]
 62
 63    # Localization.
 64    batch_shape = [2, 2]  # width (in grid points) of each state batch.
 65    # Increasing the width
 66    #  => quicker analysis (but less rel. speed-up by parallelzt., depending on NPROC)
 67    #  => worse (increased) rmse (but width 4 is only slightly worse than 1);
 68    #     if inflation is applied locally, then rmse might actually improve.
 69    localizer = nd_Id_localization(shape[::-1], batch_shape[::-1], jj, periodic=False)
 70
 71    Obs = {
 72        "M": Ny,
 73        "model": hmod,
 74        "noise": modelling.GaussRV(C=4 * np.eye(Ny)),
 75        "localizer": localizer,
 76    }
 77
 78    # Moving localization mask for smoothers:
 79    Obs["loc_shift"] = lambda ii, dt: ii  # no movement (suboptimal, but easy)
 80
 81    # Jacobian left unspecified coz it's (usually) employed by methods that
 82    # compute full cov, which in this case is too big.
 83
 84    return modelling.Operator(**Obs)
 85
 86
 87Obs = dict(time_dependent=lambda ko: obs_now(ko))
 88
 89
 90############################
 91# Other
 92############################
 93HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0, LP=LP_setup(obs_inds))
 94
 95
 96####################
 97# Suggested tuning
 98####################
 99# Reproducing Fig 7 from Sakov and Oke "DEnKF" paper from 2008.
100
101# Notes:
102# - If N<=25, then typically need to increase the dissipation
103#      to be almost sure to avoid divergence. See counillon2009.py for example.
104#    - We have not had the need to increase the dissipation parameter for the EnKF.
105# - Our experiments differ from Sakov's in the following minor details:
106#    - We use a batch width (unsure what Sakov uses).
107#    - The "EnKF-Matlab" code has a bug: it forgets to take sqrt() of the taper coeffs.
108#      This is equivalent to: R_actually_used = R_reported / sqrt(2).
109# - The boundary cells are all fixed at 0 by BCs,
110#   but are included in the state vector (amounting to 3% of the its length),
111#   and thus in RMSE calculations (which is not quite fair/optimal).
112
113# from dapper.mods.QG.sakov2008 import HMM                       # rmse.a:
114# xps += LETKF(mp=True, N=25,infl=1.04       ,loc_rad=10)        # 0.64
115# xps += LETKF(mp=True, N=25,infl='-N',xN=2.0,loc_rad=10)        # 0.66
116# xps += SL_EAKF(       N=25,infl=1.04       ,loc_rad=10)        # 0.62
117# xps += SL_EAKF(       N=25,infl=1.03       ,loc_rad=10)        # 0.58
118#
119# Iterative:
120# Yet to try: '-N' inflation, larger N, different loc_rad, and
121# smaller Lag (testing lag>3 was worse [with this loc_shift])
122# xps += iLEnKS('Sqrt',N=25,infl=1.03,loc_rad=12,nIter=3,Lag=2) # 0.59
123#
124# N = 45
125# xps += LETKF(mp=True, N=N,infl=1.02       ,loc_rad=10)        # 0.52
126# xps += LETKF(mp=True, N=N,infl='-N',xN=1.5,loc_rad=10)        # 0.51