dapper.mods.Lorenz96.spectral_obs

  1"""Settings for a "spectral" obs operator.
  2
  3Lorenz-96 is highly sensitive to large gradients.
  4Therefore, if we only observe every 4th (e.g.) state component,
  5the members might "blow up" during the forecast,
  6because the assimilation created large gradients.
  7(Of course, this will depend on R and dto).
  8Therefore, the HMM below instead uses "global obs",
  9where each observation captures information about the entire state vector.
 10The idea is that we can then remove observations, (rows of H)
 11one-by-one, to a much larger degree than for H = Identity.
 12
 13Ideally, we want the observations to be independent,
 14and possibly of the same magnitude
 15(i.e. we that the rows of H be orthonormal).
 16Furthermore, we want that each observation gives equal weight
 17("consideration") to each state component.
 18This can be shown to be equivalent to requiring that the state
 19component is equally resolved by each observation,
 20and moreover, that the magnitude (abs) of each element of H
 21be a constant (1/sqrt(Nx)).
 22
 23Can such an H be constructed/found?
 24In the 2d case: H = [1, 1; 1, -1] / sqrt(2).
 25In the 3d case: no, as can be shown by enumeration.
 26(note, however, how easy my geometric intuition was fooled.
 27Try rotating the 3-dim stensil.
 28Intuitively I thought that it would yield the 3d H of +/- 1's).
 29...
 30In fact, only in the 2^n - dimensional case is it possible
 31(our conjecture: Madleine/Patrick, based on analogy with the FFT).
 32
 33Another idea is then to evaluate the value of 40 orthogonal basis
 34functions at 40 equidistant locations
 35(corresponding to the indices of Lorenz-96).
 36This will not yield a matrix of +/- 1's,
 37but should nevertheless give nicely distributed weights.
 38
 39Note that the legendre polynomials are not orthogonal when
 40(the inner product is) evaluated on discrete, equidistant points.
 41Moreover, the actual orthogonal polynomial basis
 42(which I think goes under the name of Gram polynomials, and can be
 43constructed by qr-decomp (gram-schmidt) of a 40-dim Vandermonde matrix).
 44would in fact not be a good idea:
 45it is well-known that not only will the qr-decomp be numerically unstable,
 46the exact polynomial interpolant of 40 equidistant points
 47is subject to the "horrible" Runge phenomenon.
 48
 49Another basis is the harmonic (sine/cosine) functions. Advantages:
 50 - will be orthogonal when evaluated on 40 discrete equidistant points.
 51 - the domain of Lorenz-96 is periodic: as are sine/cosine.
 52 - (conjecture) in the 2^n dim. case, it yields a matrix of +/- 1's.
 53 - nice "spectral/frequency" interpretation of each observation.
 54Disadvatages:
 55 - 40 is not 2^n for any n
 56 - Too obvious (not very original).
 57
 58In conclusion, we will use the harmonic functions.
 59
 60
 61Update: It appears that we were wrong concerning the 2^n case.
 62That is, sine/cosine functions do not yield only +/- 1's for n>2
 63The question then remains (to be proven combinatorically?)
 64if the +/- 1 matrices exist for dim>4
 65Furthermore, experiments do not seem to indicate that I can push
 66Ny much lower than for the case H = Identity,
 67even though the rmse is a lot lower with spectral H.
 68Am I missing something?
 69"""
 70
 71import numpy as np
 72
 73import dapper.mods as modelling
 74from dapper.mods.Lorenz96.sakov2008 import Dyn, Nx, tseq
 75
 76# The (Nx-Ny) highest frequency observation modes are left out of H below.
 77# If Ny>Nx, then H no longer has independent (let alone orthogonal) columns,
 78# yet more information is gained, since the observations are noisy.
 79Ny = 12
 80
 81X0 = modelling.GaussRV(M=Nx, C=0.001)
 82
 83
 84def make_H(Ny, Nx):
 85    xx = np.linspace(-1, 1, Nx + 1)[1:]
 86    H = np.zeros((Ny, Nx))
 87    H[0] = 1 / np.sqrt(2)
 88    for k in range(-(Ny // 2), (Ny + 1) // 2):
 89        ind = 2 * abs(k) - (k < 0)
 90        H[ind] = np.sin(np.pi * k * xx + np.pi / 4)
 91    H /= np.sqrt(Nx / 2)
 92    return H
 93
 94
 95H = make_H(Ny, Nx)
 96# plt.figure(1).gca().matshow(H)
 97# plt.figure(2).gca().matshow(H @ H.T)
 98
 99Obs = {
100    "M": Ny,
101    "model": lambda x: x @ H.T,
102    "noise": modelling.GaussRV(C=0.01 * np.eye(Ny)),
103}
104
105HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0)
106
107####################
108# Obs plotting -- needs updating
109####################
110# "raw" obs plotting
111# if Ny<=Nx:
112#     Hinv = H.T
113# else:
114#     Hinv = sla.pinv2(H)
115# def yplot(y):
116#     lh = plt.plot(y @ Hinv.T,'g')[0]
117#     return lh
118
119# "implicit" (interpolated sine/cosine) obs plotting
120# Hplot = make_H(Ny,max(Ny,201))
121# Hplot_inv = Hplot.T
122# def yplot(y):
123#     x = y @ Hplot_inv.T
124#     ii = np.linspace(0,Nx-1,len(x))
125#     lh = plt.plot(ii,x,'g')[0]
126#     return lh
127
128
129####################
130# Suggested tuning
131####################
132# xps += EnKF ('Sqrt',N=40, infl=1.01)