Merge pull request #3 from UMD-AOSC/day2-dev

update nudging, hybrid, compute analysis LEs

Author: StevePny

__pycache__/class_da_system.cpython-36.pyc

@@ -73,7 +73,7 @@
 #-----------------------------------------------------------------------
 # Initialize the ensemble
 #-----------------------------------------------------------------------
-das.edim = 3 #np.int(1*xdim)
+das.edim = 2 #np.int(1*xdim)
 das.ens_bias_init = 0
 das.ens_sigma_init = 0.1
 
@@ -85,18 +85,18 @@
 I = np.identity(xdim)
 
 # Set background error covariance
-sigma_b = 1.0
+sigma_b = 1.0 #1.0
 B = I * sigma_b**2
 
 # Set observation error covariance
-sigma_r = 1.0
+sigma_r = 0.1
 R = I * sigma_r**2
 
 # Set the linear observation operator matrix as the identity by default 
 H = I
 
 # Set constant matrix for nudging
-const = 0.5
+const = 1.0
 C = I * const
 
 das.setB(B)

__pycache__/class_obs_data.cpython-36.pyc

@@ -222,7 +222,7 @@ def nudging(self,xb,yo):
     Ht = np.matrix(self.Ht)
 
     C = self.C
-    xa = xb + C*(yo - Hl*xb) 
+    xa = xb + np.dot(C,Ht)*(yo - Hl*xb) 
 
     if verbose:
       print('xb = ')
@@ -238,7 +238,7 @@ def nudging(self,xb,yo):
       print('Ht = ')
       print(Ht)
 
-    KH = Ht*C*Hl
+    KH = np.dot(np.dot(Hl,C),Ht)
 
     return xa.A1,KH
 
@@ -577,6 +577,8 @@ def HybridGain(self,Xb,yo):
 #---------------------------------------------------------------------------------------------------
 # Use a hybrid method to compute the analysis
 
+    verbose = False
+
     # Make sure inputs are matrix and column vector types
     Xb = np.matrix(Xb)
     yo = np.matrix(yo).flatten().T
@@ -590,16 +592,22 @@ def HybridGain(self,Xb,yo):
 
     # Compute ETKF analysis
     Xa_ETKF, KH_ETKF = self.ETKF(Xb,yo)
-    print('Xa_ETKF = ')
-    print(Xa_ETKF)
+    if verbose:
+      print('Xa_ETKF = ')
+      print(Xa_ETKF)
 
     # Compute 3DVar analysis
     xa_ETKF = np.mean(Xa_ETKF,axis=1)
-    print('xa_ETKF = ')
-    print(xa_ETKF)
+
+    if verbose:
+      print('xa_ETKF = ')
+      print(xa_ETKF)
+
     xa_3DVar, KH_3DVar = self._3DVar(xa_ETKF,yo)
-    print('xa_3DVar = ')
-    print(xa_3DVar)
+
+    if verbose:
+      print('xa_3DVar = ')
+      print(xa_3DVar)
 
     xa_ETKF = np.matrix(xa_ETKF).flatten().T
     xa_3DVar = np.matrix(xa_3DVar).flatten().T
@@ -611,8 +619,9 @@ def HybridGain(self,Xb,yo):
     xa_hybrid = (1-alpha)*xa_ETKF + alpha*xa_3DVar
     Xa = Xa_ETKF + np.matlib.repmat(xa_hybrid,1,edim)
 
-    print('xa_hybrid = ')
-    print(xa_hybrid)
+    if verbose:
+      print('xa_hybrid = ')
+      print(xa_hybrid)
 
     KH = (1-alpha)*KH_ETKF + alpha*KH_3DVar
 

analysis_init.py

@@ -2,7 +2,7 @@
 import pickle
 
 class obs_data:
-  def __init__(self,t=[0],pos=[0],val=[0],err=[0],bias=[],xt=[0],name='uninitialized'):
+  def __init__(self,t=[0],pos=[0],val=[0],err=[0],bias=[0],xt=[0],name='uninitialized',mu_init=[],sigma_init=[]):
     tdim = 0
     xdim = 0
     odim = 0
@@ -16,9 +16,13 @@ def __init__(self,t=[0],pos=[0],val=[0],err=[0],bias=[],xt=[0],name='uninitializ
     self.xt = np.array(xt)
     self.dep = []
     self.climvar = []
+    self.mu_init = mu_init
+    self.sigma_init = sigma_init
 
   def __str__(self):
     print(self.name)
+    print('mu_init = ', self.mu_init)
+    print('sigma_init = ', self.sigma_init)
     print('Obs Error:')
     print(self.err)
     print('Positions:')

class_da_system.py

@@ -106,10 +106,13 @@
 
   # Archive the KH matrix
   KH_history.append(deepcopy(KH))
-  KH_idx.append(i)
+  KH_idx.append(i+acyc_step)
 
 das.setKH(KH_history,KH_idx)
 
+print('Last background error covariance matrix B = ')
+print(das.getB())
+
 sv.setTrajectory(xa_history)
 sv.setName(name)
 das.setStateVector(sv)

class_obs_data.py

@@ -130,7 +130,7 @@
 
   # Archive the KH matrix
   KH_history.append(deepcopy(KH))
-  KH_idx.append(i)
+  KH_idx.append(i+acyc_step)
 
 das.setKH(KH_history,KH_idx)
 

generate_analysis_3dDet.py

@@ -0,0 +1,132 @@
+# Tutorial: "A tour of Data Assimilation methods"
+# Model: Lorenz-63
+# DA Methods: Nudging, 3D-Var, 4D-Var, Particle Filter, EnKF, Hybrid
+import numpy as np
+from class_lorenz63 import lorenz63
+from class_state_vector import state_vector
+from class_obs_data import obs_data
+from class_da_system import da_system
+from sys import argv
+from copy import deepcopy
+
+#-----------------------------------------------------------------------
+# Usage:
+#  python generate_analysis_LEs.py <method>
+#-----------------------------------------------------------------------
+
+#-----------------------------------------------------------------------
+# Read the L63 nature run
+#-----------------------------------------------------------------------
+name1 = 'x_nature'
+infile1 = name1+'_Mhist.pkl'
+sv1 = state_vector()
+sv1 = sv1.load(infile1)
+x_nature = sv1.getTrajectory()
+t_nature = sv1.getTimes()
+
+#-----------------------------------------------------------------------
+# Read the analysis
+#-----------------------------------------------------------------------
+name = 'x_analysis'
+method = argv[1]
+infile2 = name+'_'+method+'.pkl'
+das = da_system()
+das = das.load(infile2)
+sv = das.getStateVector()
+KH, khidx = das.getKH()
+print(sv)
+
+#print('KH = ')
+#print(KH)
+#print('khidx = ')
+#print(khidx)
+
+#-----------------------------------------------------------------------
+# Initialize the timesteps
+#-----------------------------------------------------------------------
+t_nature = sv1.getTimes()
+dt = (t_nature[1] - t_nature[0])
+_,xdim = np.shape(x_nature)
+acyc_step = das.acyc_step
+dtau = das.dtau
+
+#------------------------------------------------------------------
+# Propagate the TLM in time and intermittently 
+# performing a QR decomposition to rescale as necessary.
+#------------------------------------------------------------------
+rescale_interval = 2*acyc_step
+I = np.identity(xdim)
+Mhist = sv1.getMhist()
+M2hist = []
+Qhist = []
+Rhist = []
+hist_idx = []
+M2 = I
+maxit = das.maxit
+ki=0
+for i in range(maxit - acyc_step):
+  
+  # In order to find the eigenvalues of the final matrix describing
+  # the evolution of the entire trajectory, we would like to 
+  # compute the product of the M matrices to identify the
+  # stretching and contracting that results.
+  #
+  # Iterate the linear propagator (evolution matrix) 
+  # to fit the specified rescaling time interval
+
+  M2 = np.dot(Mhist[i],M2)
+
+  # If there was an analysis at this time index, then
+  # apply the contraction (I-KH) applied by the DA
+  if khidx[ki] == i:
+    M2 = np.dot(I-KH[ki],M2) 
+    ki = ki+1
+
+  # Rescale via QR decomposition
+  if (np.mod(i+1,rescale_interval)==0):
+    Q,R = np.linalg.qr(M2)
+    
+    # Store the Q and R matrices
+    M2hist.append(deepcopy(M2))
+    Qhist.append(deepcopy(Q))
+    Rhist.append(deepcopy(R))
+    hist_idx.append(i)
+    M2 = Q
+
+print('Last Q = ')
+print(Q)
+
+print('Last R = ')
+print(R)
+
+sv.rescale_interval = rescale_interval
+sv.hist_ainc = acyc_step
+sv.hist_dtau = dtau
+sv.hist_idx = hist_idx
+sv.setM2hist(M2hist)
+sv.setQhist(Qhist)
+sv.setRhist(Rhist)
+
+#-----------------------------------------------------------------------
+# Compute the Lyapunov Exponents (LEs) of the data assimilation system
+#-----------------------------------------------------------------------
+lyap_sum = np.zeros(xdim)
+rescale_cnt = 0
+for i in range(len(Rhist)):
+  R = Rhist[i]
+  rdiag = abs(np.diag(R))
+  lyap_sum = lyap_sum + np.log(rdiag) / (dtau)
+  rescale_cnt = rescale_cnt + 1
+
+lyap_avg = lyap_sum / rescale_cnt
+sv.setLEs(lyap_avg)
+
+print(sv)
+print('Lyapunov Exponents:')
+print(lyap_avg)
+
+#------------------------------------------------------------------
+# Store the nature run data
+#------------------------------------------------------------------
+outfile = name+'_LEs.pkl'
+sv.save(outfile)

generate_analysis_3dEns.py

@@ -22,8 +22,8 @@
 # Initialize the timesteps
 #-----------------------------------------------------------------------
 t_nature = sv.getTimes()
-ainc_step = 10  # (how frequently to perform an analysis)
-dtau = (t_nature[ainc_step] - t_nature[0])
+inc_step = 10  # (how frequently to perform an analysis)
+dtau = (t_nature[inc_step] - t_nature[0])
 maxit,xdim = np.shape(x_nature)
 
 #------------------------------------------------------------------
@@ -34,7 +34,7 @@
 Mhist = sv.getMhist()
 #print('Mhist = ')
 #print(Mhist)
-rescale_interval = ainc_step
+rescale_interval = inc_step
 
 M2hist = []
 Qhist = []
@@ -77,7 +77,7 @@
 print(R)
 
 sv.rescale_interval = rescale_interval
-sv.hist_ainc = ainc_step
+sv.hist_ainc = inc_step
 sv.hist_dtau = dtau
 sv.hist_idx = hist_idx
 sv.setM2hist(M2hist)

generate_analysis_LEs.py

@@ -15,7 +15,7 @@
 #--------------------------------------------------------------------------------
 # Create observation object
 #--------------------------------------------------------------------------------
-obs = obs_data(name='observe_full_state',err=[],bias=mu,val=[],pos=[])
+obs = obs_data(name='observe_full_state', mu_init=mu, sigma_init=sigma)
 
 #--------------------------------------------------------------------------------
 # Read the nature run
@@ -59,6 +59,8 @@
 obs.setDep(yo-hx)
 obs.setPos(pos)
 
+print(obs)
+
 #--------------------------------------------------------------------------------
 # Store the true and noisy observations
 #--------------------------------------------------------------------------------