-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathuntitled.txt~
100 lines (100 loc) · 3.47 KB
/
untitled.txt~
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
89
90
91
92
93
94
95
96
97
98
99
100
"def klcomp(binsize):\n",
"\n",
"\n",
" nens = 13\n",
" nEOF = 500\n",
" reg = 1\n",
"\n",
" # Load ctrl run\n",
"\n",
" out = np.load('input/CESM_ctrl_wtd_SVD.npz')\n",
" u = out['u']\n",
" s = out['s']\n",
" vt = out['vt']\n",
" lat = out['lat']\n",
" lon = out['lon']\n",
" time = out['time']\n",
" nt = out['nt']\n",
" nlat = out['nlat']\n",
" nlon = out['nlon']\n",
"\n",
" ds_TS = u.dot(np.diag(s)).dot(vt).reshape(1000,nlat,nlon)\n",
"\n",
" [nt,nlat,nlon] = np.shape(ds_TS);\n",
"\n",
" # Reshape the control run to look like a short ensemble simulation with 13 members\n",
" # New time length for these is 988 = 13*76\n",
" cnens = 13;\n",
" tdn = int(np.floor(1000./cnens)*cnens)\n",
" el = int(tdn/cnens)\n",
"\n",
" # Reshape to give an ensemble axis and transpose to make the ordering consistent\n",
" ce1 = ds_TS[:(tdn),:,:].reshape(tdn,nlat*nlon).transpose([1,0])\n",
" ce2 = ce1.reshape(nlat*nlon,el,cnens)\n",
"\n",
" # time, space, nens\n",
" ce = ce2.transpose([1,0,2])\n",
"\n",
" # Need to compute reduced-space form\n",
" [cer,uce,sce] = ens.reduce_space(ce,nEOF)\n",
"\n",
" [Cc,tCc] = ens.mk_covs(cer,np.arange(el),binsize)\n",
" [tdc,_,_] = Cc.shape\n",
" Cvc = Cc.reshape(tdc,nEOF**2).T\n",
"\n",
" # No smoothing for mean\n",
" [Ccm,tCcm] = ens.mk_covs(cer,np.arange(el),1)\n",
" [tdcm,_,_] = Ccm.shape\n",
" Cvcm = Ccm.reshape(tdcm,nEOF**2).T\n",
" Cmc = Cvcm.mean(axis=1, keepdims=True).reshape(nEOF,nEOF)\n",
"\n",
" m0 = np.zeros(nEOF)\n",
" m1 = np.zeros(nEOF)\n",
"\n",
" kldc = np.empty(tdc)\n",
" for ii in np.arange(tdc):\n",
" kldc[ii] = ens.KLdiv(Cc[ii,:,:]+reg*np.eye(nEOF),Cmc+reg*np.eye(nEOF),m0,m1)\n",
"\n",
" \n",
" \n",
" \n",
"\n",
" # Now for LME\n",
" \n",
" out = np.load('input/CESM_LME_all13_wtd_SVD.npz')\n",
" u = out['u']\n",
" s = out['s']\n",
" vt = out['vt']\n",
" lat = out['lat']\n",
" lon = out['lon']\n",
" time = out['time']\n",
" nt = out['nt']\n",
" nlat = out['nlat']\n",
" nlon = out['nlon']\n",
" nens = out['nens']\n",
"\n",
" # EOF x time*nens\n",
" datr = (vt[:nEOF,:]*s[:nEOF,None])\n",
"\n",
" # reshaped into timexEOFxnens\n",
" datrr = datr.reshape(nEOF,nt,nens).transpose(1,0,2)\n",
"\n",
" # Get time-varying reduced-space covariances\n",
" [C,tC] = ens.mk_covs(datrr,time,binsize)\n",
" [td,_,_] = C.shape\n",
" Cv = C.reshape(td,nEOF**2).T\n",
"\n",
" # No smoothing for mean\n",
" [Cm,tCm] = ens.mk_covs(cer,np.arange(el),1)\n",
" [tdm,_,_] = Ccm.shape\n",
" Cvm = Ccm.reshape(tdcm,nEOF**2).T\n",
" Cm = Cvm.mean(axis=1, keepdims=True).reshape(nEOF,nEOF)\n",
"\n",
" m0 = np.zeros(nEOF)\n",
" m1 = np.zeros(nEOF)\n",
"\n",
" kld = np.empty(td)\n",
" for ii in np.arange(td):\n",
" kld[ii] = ens.KLdiv(C[ii,:,:]+reg*np.eye(nEOF),Cmc+reg*np.eye(nEOF),m0,m1)\n",
"\n",
" return tC, kld, tCc, kldc"