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

Day2 dev

Author: StevePny

README.md

@@ -2,3 +2,15 @@
 
 This is a 'hands-on' tutorial for the RIKEN International School on Data Assimilation (RISDA2018).
 http://www.data-assimilation.riken.jp/risda2018/
+
+This software package is based on routines developed for the publication:
+
+Mathematical foundations of hybrid data assimilation from a synchronization perspective
+Penny, S.G., Chaos 27, 126801 (2017); https://doi.org/10.1063/1.5001819
+
+http://aip.scitation.org/doi/full/10.1063/1.5001819
+
+Instructions:
+This package uses the anaconda python distribution and for plotting uses plotly.
+
+

__pycache__/class_da_system.cpython-36.pyc

@@ -3,6 +3,21 @@
 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
+
+#-----------------------------------------------------------------------
+# Usage:
+# python analysis_init.py {method}
+#
+# where {method} is any one of:
+#  skip
+#  nudging
+#  OI
+#  3DVar
+#  ETKF
+#  PF
+#  Hybrid
+#-----------------------------------------------------------------------
 
 #-----------------------------------------------------------------------
 # Read the L63 nature run
@@ -23,12 +38,15 @@
 #-----------------------------------------------------------------------
 # Try reducing the observed dimensions
 #-----------------------------------------------------------------------
-#obs.reduceDim([0])  # x only
-#obs.reduceDim([1])  # y only
-#obs.reduceDim([2])  # z only
-#obs.reduceDim([0,1])  # x and y only
-#obs.reduceDim([1,2])  # y and z only
-#obs.reduceDim([0,2])  # z and x only
+yp = [0,1,2]
+#yp = [0]    # x only
+#yp = [1]    # y only
+#yp = [2]    # z only
+#yp = [0,1]  # x and y only
+#yp = [1,2]  # y and z only
+#yp = [0,2]  # z and x only
+if len(yp) < xdim:
+  obs.reduceDim(yp)
 
 y_obs = obs.getVal()
 y_pts = obs.getPos()
@@ -52,6 +70,13 @@
 das.t = sv.getTimes()
 das.t0 = das.t[0]
 
+#-----------------------------------------------------------------------
+# Initialize the ensemble
+#-----------------------------------------------------------------------
+das.edim = 3 #np.int(1*xdim)
+das.ens_bias_init = 0
+das.ens_sigma_init = 0.1
+
 #-----------------------------------------------------------------------
 # Initialize the error covariances B and R, and the linearized 
 # observation operator H
@@ -60,7 +85,7 @@
 I = np.identity(xdim)
 
 # Set background error covariance
-sigma_b = 0.9
+sigma_b = 1.0
 B = I * sigma_b**2
 
 # Set observation error covariance
@@ -70,16 +95,26 @@
 # Set the linear observation operator matrix as the identity by default 
 H = I
 
-print('B = ')
-print(B)
-print('R = ')
-print(R)
-print('H = ')
-print(H)
+# Set constant matrix for nudging
+const = 0.5
+C = I * const
 
 das.setB(B)
 das.setR(R)
 das.setH(H)
+das.setC(C)
+
+# Update the matrices to fit the reduced observation dimension
+if len(yp) < xdim:
+  das.reduceYdim(yp)
+
+print('B = ')
+print(das.getB())
+print('R = ')
+print(das.getR())
+print('H = ')
+print(das.getH())
+
 
 #-----------------------------------------------------------------------
 # Initialize the timesteps
@@ -104,10 +139,12 @@
 # Choose DA method:
 #-----------------------------------------------------------------------
 
+method = argv[1]
+
 #-----------
 # Test basic functionality
 #-----------
-method='skip'
+#method='skip'
 
 #-----------
 # 3D methods

__pycache__/class_lorenz63.cpython-36.pyc

@@ -29,6 +29,12 @@ def __init__(self,x0=[],yo=[],t0=0,dt=0,alpha=0.5,state_vector=[],obs_data=[],ac
     self.SqrtB = []
     self.state_vector = state_vector
     self.obs_data = obs_data
+    self.method = ''
+    self.KH = []
+    self.khidx = []
+    self.edim = 0
+    self.das_bias_init = 0
+    self.das_sigma_init = 1
 
   def __str__(self):
     print('xdim = ', self.xdim)
@@ -51,6 +57,8 @@ def __str__(self):
     print(self.state_vector)
     print('obs_data = ')
     print(obs_data)
+    print('method = ')
+    print(self.method)
     return 'type::da_system'
 
   def setMethod(self,method):
@@ -80,6 +88,12 @@ def update(self,B=[0],R=[0],H=[0],t=[0],x0=[0]):
     self.t = t
     self.x0 = x0
 
+  def getC(self):
+    return self.C
+
+  def setC(self,C):
+    self.C = np.matrix(C)
+
   def getB(self):
     return self.B
 
@@ -119,12 +133,22 @@ def setH(self,H):
 #   if (nc != self.xdim):
 #     error('H must be ydim x xdim, but instead H is %d x %d'%(nr,nc))
 
+  def getKH(self):
+    return self.KH, self.khidx
+
+  def setKH(self,KH,khidx):
+    self.KH = KH
+    self.khidx = khidx
+
   def reduceYdim(self,yp):
 #   print('reduceYdim:')
 #   print('yp = ', yp)
     self.ydim = len(yp)
     self.setH(self.H[yp,:])
-    self.setR(self.R[yp,yp])
+    R = self.R
+    R = R[yp,:]
+    R = R[:,yp]
+    self.setR(R)
 
   def compute_analysis(self,xb,yo,params=[0]):
     # (params needed for 4D-Var)
@@ -190,16 +214,30 @@ def nudging(self,xb,yo):
 #---------------------------------------------------------------------------------------------------
 # Use observations at predefined points to drive the model system to the observed nature system
 
+    verbose = False
+
     xb = np.matrix(xb).flatten().T
     yo = np.matrix(yo).flatten().T
+    Hl = np.matrix(self.H)
+    Ht = np.matrix(self.Ht)
 
-    const = np.diagonal(self.B)
+    C = self.C
+    xa = xb + C*(yo - Hl*xb) 
 
-    xa = xb + const*(yo - xb) 
+    if verbose:
+      print('xb = ')
+      print(xb)
+      print('C = ')
+      print(C)
+      print('yo = ')
+      print(yo)
+      print('Hl = ')
+      print(Hl)
+      print('xa = ')
+      print(xa)
+      print('Ht = ')
+      print(Ht)
 
-    C = np.diag(const)
-    Hl = self.H
-    Ht = self.Ht
     KH = Ht*C*Hl
 
     return xa.A1,KH
@@ -255,10 +293,10 @@ def _3DVar(self,xb,yo):
 
     # 'preconditioning with B'
     I = np.identity(xdim)
-    BHt = B*Ht
-    BHtRinv = BHt*Rinv
-    A = I + BHtRinv*Hl
-    b1 = xb + BHtRinv*yo
+    BHt = np.dot(B,Ht)
+    BHtRinv = np.dot(BHt,Rinv)
+    A = I + np.dot(BHtRinv,Hl)
+    b1 = xb + np.dot(BHtRinv,yo)
 
     # Use minimization algorithm to minimize cost function:
     xa,ierr = sp.sparse.linalg.cg(A,b1,x0=xb,tol=1e-05,maxiter=1000) 

__pycache__/class_obs_data.cpython-36.pyc

@@ -51,8 +51,8 @@ def setTrajectory(self,states):
     x_avg = np.zeros(self.xdim)
     x_std = np.zeros(self.xdim)
     for i in range(self.xdim):
-      x_avg[i] = np.mean(states[:,i])
-      x_std[i] = np.std(states[:,i])
+      x_avg[i] = np.nanmean(states[:,i])
+      x_std[i] = np.nanstd(states[:,i])
     self.clim_mean = x_avg
     self.clim_std = x_std
 

__pycache__/class_state_vector.cpython-36.pyc

@@ -0,0 +1,2 @@
+rm *_*.pkl
+rm *.html

analysis_init.py

@@ -57,33 +57,17 @@
 #-----------------------------------------------------------------------
 # Choose DA method:
 #-----------------------------------------------------------------------
-
-method = das.getMethod()  # (use default)
-
-#-----------
-# Test basic functionality
-#-----------
-#method='skip'
-
-#-----------
-# 3D methods
-#-----------
-# Nudging
-#method='nudging'
-# OI
-#method='OI'
-# 3D-Var
-#method='3DVar'
-
-#das.setMethod(method)
+method = das.getMethod()
 
 #-----------------------------------------------------------------------
 # Conduct data assimilation process
 #-----------------------------------------------------------------------
 #
-xa = sv.x0
+xa = das.x0
 xa_history = np.zeros_like(x_nature)
+xa_history[:] = np.nan
 KH_history = []
+KH_idx = []
 for i in range(0,maxit-acyc_step,acyc_step):
 
   #----------------------------------------------
@@ -104,9 +88,6 @@
   yo = y_obs[i+acyc_step,:]
   yp = y_pts[i+acyc_step,:]
 
-# if (len(yp) < xdim): 
-#   das.reduceYdim(yp)
- 
   #----------------------------------------------
   # Compute analysis
   #----------------------------------------------
@@ -116,22 +97,18 @@
 # print('KH = ')
 # print(KH)
 
-  # Archive the analysis
-  xa_history[i+acyc_step,:] = xa
   # Fill in the missing timesteps with the forecast from the previous analysis IC's
   xa_history[i:i+acyc_step,:] = xf_4d[0:acyc_step,:]
+  # Archive the analysis
+  xa_history[i+acyc_step,:] = xa
 
 # print('xa_history[i:i+acyc_step+1,:] = ', xa_history[i:i+acyc_step+1,:])
 
   # Archive the KH matrix
   KH_history.append(deepcopy(KH))
+  KH_idx.append(i)
 
-#--------------------------------------------------------------------------------
-# Fill in unobserved dimensions (for plotting)
-#--------------------------------------------------------------------------------
-#fillValue=0.0
-#obs.fillDim(fillValue)
-#das.setObsData(obs)
+das.setKH(KH_history,KH_idx)
 
 sv.setTrajectory(xa_history)
 sv.setName(name)