-
Notifications
You must be signed in to change notification settings - Fork 0
/
test_soar.py
executable file
·88 lines (79 loc) · 2.44 KB
/
test_soar.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#!/usr/bin/env python
import numpy as np
import misc
import data_assimilation as DA
import config as p
import matplotlib.pyplot as plt
from scipy.fftpack import fft
plt.switch_backend('Agg')
plt.figure(figsize=(8, 8))
nx = p.nx
v = misc.fourier_basis(nx)
###SOAR process with parameters f1 and f2
####x(t) = f1 x(t-1) + f2 x(t-2) + e(t)
nsample = 10000
x = np.zeros((nx, nsample))
f1 = 0.6
f2 = -0.0
sig = 1
for n in range(nsample):
x[0, n] = np.random.normal(0, sig)
x[1, n] = f1*x[0, n] + np.random.normal(0, sig)
x[nx-1, n] = f1*x[0, n] + np.random.normal(0, sig)
# for i in range(int(nx/2)-2):
for i in range(nx-2):
x[i+2, n] = f1*x[i+1, n] + f2*x[i, n] + np.random.normal(0, sig)
# for i in range(nx-1, int(nx/2)+1, -1):
# x[i-2, n] = f1*x[i-1, n] + f2*x[i, n] + np.random.normal(0, sig)
R = misc.error_covariance(x)
L = -1/np.log(f1)
sig1 = np.sqrt(sig**2/(1-f1**2))
R1 = DA.R_matrix(nx, 1, p.obs_ind, sig1, L, 0, 1)
print(np.mean(np.diag(R)))
print(np.mean(np.diag(R1)))
ax = plt.subplot(221)
##spectrum from simulated R
wo = np.diag(np.dot(v.T, np.dot(R, v)))
ax.plot(np.sqrt(wo[::2]), 'r')
##spectrum from theory
ii = 2*np.pi*np.mgrid[0:int(nx/2)]/nx
# wo2 = sig**2/(1 + f1**2 - 2*f1*np.cos(ii)) ##AR1
wo2 = sig**2/(1 + f1**2 + f2**2 - 2*f1*(1-f2)*np.cos(ii) - 2*f2*np.cos(2*ii)) ##AR2
ax.plot(np.sqrt(wo2), 'b')
##spectrum from true R
wo1 = np.diag(np.dot(v.T, np.dot(R1, v)))
ax.plot(np.sqrt(wo1[::2]), 'k')
ax.set_ylim(0, 3)
ax.set_xlabel('wavenumber')
ax = plt.subplot(222)
ax.plot(R[20, :], 'r') ##simulated R
jj = np.mgrid[0:nx] - 20
R2 = sig1**2 * np.exp(-np.abs(jj)/L)
ax.plot(R2, 'b') ##from theory
ax.plot(R1[20, :], 'k') ##true R
ax.set_xlim(0, nx)
ax.set_ylim(-1, 2)
ax.set_xlabel('distance')
ax = plt.subplot(223)
ax.contourf(R, np.arange(-2, 2.1, 0.1), cmap='seismic')
ax = plt.subplot(224)
ax.contourf(R1, np.arange(-2, 2.1, 0.1), cmap='seismic')
###spectral distribution
# ax = plt.subplot(223)
# for L in [0, 2, 5, 10]:
# R = DA.R_matrix(nx, 1, p.obs_ind, 1, L, 0, 1)
# wo = np.diag(np.dot(v.T, np.dot(R, v)))
# ax.plot(np.sqrt(wo[::2]), label='L = {}'.format(L))
# ax.legend(fontsize=10)
# ax.set_ylim(0, 3)
# ax.set_xlabel('wavenumber')
###correlation length scale
# ax = plt.subplot(224)
# for L in [0, 2, 5, 10]:
# R = DA.R_matrix(nx, 1, p.obs_ind, 1, L, 0, 1)
# ax.plot(R[0, :], label='L = {}'.format(L))
# ax.legend(fontsize=10)
# ax.set_xlim(0, 20)
# ax.set_ylim(-1, 2)
# ax.set_xlabel('distance')
plt.savefig('1.pdf')