in flight changes

class_da_system.py

@@ -7,7 +7,7 @@
 
 class da_system:
 
-  def __init__(self,x0=[0],yo=[0],t0=0,dt=0,alpha=0.5,state_vector=[0],obs_data=[0]):
+  def __init__(self,x0=[],yo=[],t0=0,dt=0,alpha=0.5,state_vector=[],obs_data=[],acyc_step=10):
     self.xdim = np.size(x0)
     self.ydim = np.size(yo)
     self.edim = 1
@@ -16,15 +16,19 @@ def __init__(self,x0=[0],yo=[0],t0=0,dt=0,alpha=0.5,state_vector=[0],obs_data=[0
     self.dt = dt
     self.X0 = x0
     self.t = t0
-    self.ainc = 1
+    self.acyc_step = acyc_step
+    self.dtau = dt*acyc_step
+    self.fcst_step = acyc_step
+    self.fcst_dt = dt
+    self.maxit = 0
     self.B = np.matrix(np.identity(self.xdim))
     self.R = np.matrix(np.identity(self.ydim))
     self.H = np.matrix(np.identity(self.xdim))
     self.Ht = (self.H).transpose()
     self.alpha = alpha
-    self.SqrtB = sp.linalg.sqrtm(self.B)
-    self.obs_data = obs_data
+    self.SqrtB = []
     self.state_vector = state_vector
+    self.obs_data = obs_data
 
   def __str__(self):
     print('xdim = ', self.xdim)
@@ -33,10 +37,20 @@ def __str__(self):
     print('t0 = ', self.t0)
     print('dt = ', self.dt)
     print('t  = ', self.t)
+    print('acyc_step = ', self.acyc_step)
+    print('dtau = ', self.dtau)
+    print('fcst_step = ', self.fcst_step)
+    print('fcst_dt  = ', self.fcst_dt)
     print('B = ')
     print(self.B)
     print('R = ')
     print(self.R)
+    print('H = ')
+    print(self.H)
+    print('state_vector = ')
+    print(self.state_vector)
+    print('obs_data = ')
+    print(obs_data)
     return 'type::da_system'
 
   def setMethod(self,method):
@@ -113,42 +127,56 @@ def reduceYdim(self,yp):
     self.setR(self.R[yp,yp])
 
   def compute_analysis(self,xb,yo,params=[0]):
+    # (params needed for 4D-Var)
     method = self.method
     if method == 'skip':
       xa = xb 
+      KH = np.identity(self.xdim)
     elif method == 'nudging':
-      xa = self.nudging(xb,yo)
+      xa,KH = self.nudging(xb,yo)
     elif method == 'OI':
-      xa = self.OI(xb,yo)
+      xa,KH = self.OI(xb,yo)
     elif method == '3DVar' or method == '3D-Var':
-      xa = self._3DVar(xb,yo) 
+      xa,KH = self._3DVar(xb,yo) 
     elif method == 'ETKF' or method == 'EnKF':
-      xa = self.ETKF(xb,yo)
+      xa,KH = self.ETKF(xb,yo)
     elif method == 'PF':
-      xa = self.PF(xb,yo)
+      xa,KH = self.PF(xb,yo)
     elif method == 'Hybrid':
-      xa = self.HybridGain(xb,yo) 
+      xa,KH = self.HybridGain(xb,yo) 
 #   elif method == '4DVar' or method == '4D-Var':
-#     xa = self._4DVar(xb,yo)
+#     xa,KH = self._4DVar(xb,yo)
 #   elif method == '4DEnVar':
-#     xa = self._4DEnVar(xb,yo)
+#     xa,KH = self._4DEnVar(xb,yo)
 #   elif method == '4DETKF':
-#     xa = self._4DETKF(xb,yo) 
+#     xa,KH = self._4DETKF(xb,yo) 
     else:
       print('compute_analysis:: Unrecognized DA method.')
       raise SystemExit
-    return xa
+    return xa,KH
 
   def initEns(self,x0,mu=0,sigma=0.1,edim=4):
     xdim = len(x0)
     x0 = np.matrix(x0).flatten().T
-    xrand = np.random.normal(mu,sigma,(xdim,edim))
+    mu = np.matrix(mu).flatten().T
+    Xrand = np.random.normal(0,sigma,(xdim,edim))
+    Xrand = np.matrix(Xrand)
+#   print('Xrand = ')
+#   print(Xrand)
+    # Remove mean to make sure it is properly centered at 0
+    # (add bias if included)
+    rand_mean = np.mean(Xrand,axis=1) + mu
+#   print('rand_mean = ')
+#   print(rand_mean)
+    rmat = np.matlib.repmat(rand_mean,1,edim)
+    Xrand = Xrand - rmat
 #   print('xrand = ')
 #   print(xrand)
+    # add perturbations to x0
     rmat = np.matlib.repmat(x0, 1, edim)
 #   print('rmat = ')
 #   print(rmat)
-    X0 = np.matrix(rmat + xrand)
+    X0 = np.matrix(rmat + Xrand)
     return X0
 
 # def init4D(self):

class_lorenz63.py

@@ -3,7 +3,7 @@
 import plotly.offline as py
 import plotly.graph_objs as go
 from plotly import tools
-#from copy import deepcopy
+from copy import deepcopy
 
 #------------------------------------------------------------------
 # Define functions
@@ -14,9 +14,10 @@ def f(state, t, sigma, rho, beta):
 
 def Ja(state, t, sigma, rho, beta):
   x, y, z = state  # unpack the state vector
-  J = [[-sigma, sigma,   0   ]
-       [ rho-z,    -1,  -x   ]
+  J = [[-sigma, sigma,   0   ],
+       [ rho-z,    -1,  -x   ],
        [     y,     x,  -beta]]  
+  J = np.matrix(J)
   return J
 
 def Jfd(state0,state1,params):
@@ -125,6 +126,44 @@ def run(self, state0, t):
     states = odeint(f, state0, t, args=(self.sigma, self.rho, self.beta))
     return states
 
+  #------------------------------------------------------------------
+  # Compute approximate TLM with I + Df(x0)*dt
+  #------------------------------------------------------------------
+  def compute_TLMa(self, states, t):
+
+    print('states = ')
+    print(states)
+
+    print('times = ')
+    print(t)
+
+    nr,nc = np.shape(states)
+    I = np.identity(nc)
+
+    # Compute Jacobian / linear propagator for each timestep
+    sigma,rho,beta = self.params
+    maxit = len(t)
+    Mhist=[]
+    for i in range(maxit):
+      if (i < maxit-1):
+        dt = t[i+1] - t[i]
+      else:
+        dt = t[-1] - t[-2]
+
+      # Evaluate Jacobian
+      Df = Ja(states[i,:], t[i], sigma, rho, beta)
+
+#     print('Df = ')
+#     print(Df)
+
+      # Compute approximate linear propagator
+      # (Note this is a poor approximation of the matrix integral, 
+      # we would prefer a geometric integrator, e.g. the Magnus expansion)
+      M = I + Df*dt  
+      Mhist.append(deepcopy(M))
+
+    return Mhist
+
   #------------------------------------------------------------------
   # Compute approximate Jacobian with finite differences
   #------------------------------------------------------------------
@@ -188,10 +227,10 @@ def plot(self, states,cvec,outfile='l63-3d',plot_title='Lorenz 63 attractor'):
     data = [trace]
 
     layout = dict(
-      width=800,
+      width=1200,
       height=700,
       autosize=False,
-      title=self.title,
+      title=plot_title,
       scene=dict(
           xaxis=dict(
               gridcolor='rgb(255, 255, 255)',

class_obs_data.py

@@ -23,8 +23,8 @@ def __str__(self):
     print(self.err)
     print('Positions:')
     print(self.pos)
-    print('Time interval:')
-    print(self.t)
+#   print('Time interval:')
+#   print(self.t)
     print('Values:')
     print(self.val)
     return self.name

class_state_vector.py

@@ -16,6 +16,8 @@ def __init__(self,params=[0],x0=[0],t=[0],name='uninitialized'):
     self.Mhist = []
     self.Qhist = []
     self.Rhist = []
+    self.hist_ainc = 0
+    self.hist_dtau = 0
     self.hist_idx = []
     self.rescale_interval = 1
     self.TLM = []
@@ -36,6 +38,9 @@ def __str__(self):
     print(self.clim_std)
     return self.name
 
+  def setName(self,name):
+    self.name = name
+
   def getTrajectory(self):
     return self.trajectory
 
@@ -75,6 +80,12 @@ def getMhist(self):
   def setMhist(self,Mhist):
     self.Mhist = Mhist
 
+  def getM2hist(self):
+    return self.M2hist
+
+  def setM2hist(self,M2hist):
+    self.M2hist = M2hist
+
   def getQhist(self):
     return self.Qhist
 

generate_analysis_3dDet.py

@@ -6,30 +6,48 @@
 from class_state_vector import state_vector
 from class_obs_data import obs_data
 from class_da_system import da_system
+from copy import deepcopy
 
 #-----------------------------------------------------------------------
 # Read the da system object
 #-----------------------------------------------------------------------
-name = 'x_analysis_init'
-infile = name+'.pkl'
+name = 'x_analysis'
+infile = name+'_init.pkl'
+das = da_system()
+das = das.load(infile)
+
+print(das)
+
+#-----------------------------------------------------------------------
+# Get the nature run trajectory
+#-----------------------------------------------------------------------
 sv = das.getStateVector()
 x_nature = sv.getTrajectory()
 
+#-----------------------------------------------------------------------
+# Get the L63 observations via the obs_data object
+#-----------------------------------------------------------------------
+obs = das.getObsData()
+y_obs = obs.getVal()
+y_pts = obs.getPos()
+y_err = obs.getErr()
+print('y_obs = ')
+print(y_obs[0,:])
+print('y_pts = ')
+print(y_pts[0,:])
+
 #-----------------------------------------------------------------------
 # Initialize the timesteps
 #-----------------------------------------------------------------------
 t_nature = sv.getTimes()
-ainc_step = das.ainc  # (how frequently to perform an analysis)
+acyc_step = das.acyc_step  # (how frequently to perform an analysis)
 dtau = das.dtau
-tsteps= das.tsteps
 dt = das.dt
+fcst_step= das.fcst_step
+fcst_dt = das.fcst_dt
 maxit = das.maxit
 xdim = das.xdim
-
-#-----------------------------------------------------------------------
-# Get the L63 observations via the obs_data object
-#-----------------------------------------------------------------------
-obs = das.getObsData()
+ydim = das.ydim
 
 #-----------------------------------------------------------------------
 # Initialize the model
@@ -43,23 +61,21 @@
 method = das.getMethod()  # (use default)
 
 #-----------
-# Session 0:
 # Test basic functionality
 #-----------
 #method='skip'
 
 #-----------
-# Session 2:
 # 3D methods
 #-----------
 # Nudging
 #method='nudging'
 # OI
-method='OI'
+#method='OI'
 # 3D-Var
 #method='3DVar'
 
-das.setMethod(method)
+#das.setMethod(method)
 
 #-----------------------------------------------------------------------
 # Conduct data assimilation process
@@ -68,31 +84,28 @@
 xa = sv.x0
 xa_history = np.zeros_like(x_nature)
 KH_history = []
-for i in range(0,maxit,ainc_step):
+for i in range(0,maxit-acyc_step,acyc_step):
 
   #----------------------------------------------
   # Run forecast model for this analysis cycle:
   #----------------------------------------------
-  t = np.arange(t_nature[i],t_nature[i]+dtau+dt,dt)
+  t = np.arange(t_nature[i],t_nature[i+acyc_step]+dt,dt)
 # print('t = ', t)
+# print('t_nature[i+acyc_step] = ', t_nature[i+acyc_step])
+
   # Run the model
-  xf_4D =  l63.run(xa,t) 
+  xf_4d =  l63.run(xa,t) 
   # Get last timestep of the forecast
-  xf = xf_4D[-1,:] 
+  xf = xf_4d[-1,:] 
 
   #----------------------------------------------
   # Get the observations for this analysis cycle
   #----------------------------------------------
-  yo = y_obs[i,:]
-  yp = y_pts[i,:]
+  yo = y_obs[i+acyc_step,:]
+  yp = y_pts[i+acyc_step,:]
 
-  #----------------------------------------------
-  # Update the error covariances
-  #----------------------------------------------
-  das.setB(sigma_b**2*I)
-  das.setR(sigma_r**2*I)
-  das.setH(I)
-  das.reduceYdim(yp)
+# if (len(yp) < xdim): 
+#   das.reduceYdim(yp)
 
   #----------------------------------------------
   # Compute analysis
@@ -104,10 +117,14 @@
 # print(KH)
 
   # Archive the analysis
-  xa_history[i,:] = xa
+  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,:]
+
+# print('xa_history[i:i+acyc_step+1,:] = ', xa_history[i:i+acyc_step+1,:])
 
   # Archive the KH matrix
-  KH_history.append(KH)
+  KH_history.append(deepcopy(KH))
 
 #--------------------------------------------------------------------------------
 # Fill in unobserved dimensions (for plotting)
@@ -117,7 +134,10 @@
 #das.setObsData(obs)
 
 sv.setTrajectory(xa_history)
+sv.setName(name)
 das.setStateVector(sv)
 
-outfile='x_analysis_'+method+'.pkl'
+outfile=name+'_'+method+'.pkl'
 das.save(outfile)
+
+print(das)