From 8af51c033c2274a90c964f34d96a1be6e153a38d Mon Sep 17 00:00:00 2001 From: System Administrator Date: Mon, 22 Jan 2018 10:45:50 +0900 Subject: [PATCH] updates for configuration --- class_da_system.py | 144 ++++++++++++++++++++++++++++------ class_lorenz63.py | 12 +-- generate_analysis_3dEns.py | 28 ++++++- generate_analysis_4dDet.py | 1 + generate_analysis_4dEns.py | 1 + generate_analysis_init.py | 155 +++++++++++++++++++++++++++++++++++++ generate_nature_run.py | 2 +- plot_analysis_vs_nature.py | 2 +- runall_tutorial_1.sh | 5 ++ 9 files changed, 318 insertions(+), 32 deletions(-) create mode 100644 generate_analysis_init.py create mode 100644 runall_tutorial_1.sh diff --git a/class_da_system.py b/class_da_system.py index de07fb7..4d63139 100755 --- a/class_da_system.py +++ b/class_da_system.py @@ -165,7 +165,7 @@ def initEns(self,x0,mu=0,sigma=0.1,edim=4): # 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 + rand_mean = np.mean(Xrand,axis=1) - mu # print('rand_mean = ') # print(rand_mean) rmat = np.matlib.repmat(rand_mean,1,edim) @@ -176,7 +176,7 @@ def initEns(self,x0,mu=0,sigma=0.1,edim=4): rmat = np.matlib.repmat(x0, 1, edim) # print('rmat = ') # print(rmat) - X0 = np.matrix(rmat + Xrand) + X0 = np.matrix(Xrand + rmat) return X0 # def init4D(self): @@ -276,6 +276,7 @@ def _3DVar(self,xb,yo): def ETKF(self,Xb,yo): #--------------------------------------------------------------------------------------------------- # Use ensemble of states to estimate background error covariance + verbose=False # Make sure inputs are matrix and column vector types Xb = np.matrix(Xb) @@ -291,12 +292,12 @@ def ETKF(self,Xb,yo): Hl = self.H Yb = np.matrix(np.zeros([ydim,edim])) for i in range(edim): - Yb[:,i] = Hl*Xb[:,i] + Yb[:,i] = np.dot(Hl,Xb[:,i]) # Convert ensemble members to perturbations xm = np.mean(Xb,axis=1) - Xb = Xb - np.matlib.repmat(xm, 1, edim) ym = np.mean(Yb,axis=1) + Xb = Xb - np.matlib.repmat(xm, 1, edim) Yb = Yb - np.matlib.repmat(ym, 1, edim) # Compute R^{-1} @@ -304,34 +305,127 @@ def ETKF(self,Xb,yo): Rinv = np.linalg.inv(R) # Compute the weights - Ybt = Yb.T - C = Ybt*Rinv + #---- + # stage(4) Now do computations for lxk Yb matrix + # Compute Yb^T*R^(-1) + #---- + Ybt = np.transpose(Yb) + C = np.dot(Ybt,Rinv) + + if verbose: + print ('C = ') + print (C) + + #---- + # stage(5) Compute eigenvalue decomposition for Pa + # Pa = [(k-1)I/rho + C*Yb]^(-1) + #---- I = np.identity(edim) - rho = 1.0 - eigArg = (edim-1)*I/rho + C*Yb + rho = 1.05 #1.0 + eigArg = (edim-1)*I/rho + np.dot(C,Yb) + lamda,P = np.linalg.eigh(eigArg) + + if verbose: + print ('lamda = ') + print (lamda) + print ('P = ') + print (P) + Linv = np.diag(1.0/lamda) - PLinv = P*Linv - Pt = P.T - Pa = PLinv*Pt + if verbose: + print ('Linv = ') + print (Linv) + + PLinv = np.dot(P,Linv) + + if verbose: + print ('PLinv = ') + print (PLinv) + + Pt = np.transpose(P) + + if verbose: + print ('Pt = ') + print (Pt) + + Pa = np.dot(PLinv, Pt) + + if verbose: + print ('Pa = ') + print (Pa) + + #---- + # stage(6) Compute matrix square root + # Wa = [(k-1)Pa]1/2 + #---- Linvsqrt = np.diag(1/np.sqrt(lamda)) - PLinvsqrt = P*Linvsqrt - Wa = np.sqrt(edim-1) * PLinvsqrt*Pt - d = yo - ym - wm = Pa*(C*d) + if verbose: + print ('Linvsqrt = ') + print (Linvsqrt) + + PLinvsqrt = np.dot(P,Linvsqrt) + + if verbose: + print ('PLinvsqrt = ') + print (PLinvsqrt) + + Wa = np.sqrt((edim-1)) * np.dot(PLinvsqrt,Pt) + + if verbose: + print ('Wa = ') + print (Wa) + + #---- + # stage(7) Transform back + # Compute the mean update + # Compute wabar = Pa*C*(yo-ybbar) and add it to each column of Wa + #---- + d = yo-ym + Cd = np.dot(C,d) + + if verbose: + print ('Cd = ') + print (Cd) + + wm = np.dot(Pa,Cd) #k x 1 + + if verbose: + print ('wm = ') + print (wm) + + # Add the same mean vector wm to each column +# Wa = Wa + wm[:,np.newaxis] #STEVE: make use of python broadcasting to add same vector to each column Wa = Wa + np.matlib.repmat(wm, 1, edim) + if verbose: + print ('Wa = ') + print (Wa) + + #---- + # stage(8) + # Compute the perturbation update + # Multiply Xb (perturbations) by each wa(i) and add xbbar + #---- + # Add the same mean vector wm to each column - Xa = Xb*Wa + np.matlib.repmat(xm, 1, edim) +# Xa = np.dot(Xb,Wa) + xm[:,np.newaxis] + Xa = np.dot(Xb,Wa) + np.matlib.repmat(xm, 1, edim) + + if verbose: + print ('Xa = ') + print (Xa) # Compute KH: - IpYbtRinvYb = ((edim-1)/rho)*I + Ybt*Rinv*Yb - IpYbtRinvYb_inv = IpYbtRinvYb.I - K = Xb*IpYbtRinvYb_inv*Ybt*Rinv - KH = K*Hl + RinvYb = np.dot(Rinv,Yb) + IpYbtRinvYb = ((edim-1)/rho)*I + np.dot(Ybt,RinvYb) + IpYbtRinvYb_inv = np.linalg.inv(IpYbtRinvYb) + YbtRinv = np.dot(Ybt,Rinv) + K = np.dot( Xb, np.dot(IpYbtRinvYb_inv,YbtRinv) ) + KH = np.dot(K,Hl) return Xa, KH @@ -350,19 +444,23 @@ def ETKF(self,Xb,yo): # R_4d = self.R ! may need list of R matrices if observations are non-uniform over time # H_4d = self.H ! may need list of H matrices if observations are non-uniform over time # M_4d = self.M ! input list of TLMs for each timestep +# +# (work in progress) #--------------------------------------------------------------------------------------------------- # def _4DEnVar(self,Xb_4d,yo_4d): #--------------------------------------------------------------------------------------------------- # Use ensemble of states over a time window to estimate temporal correlations - +# +# (work in progress) #--------------------------------------------------------------------------------------------------- # def _4DETKF(self,Xb_4d,yo_4d): #--------------------------------------------------------------------------------------------------- # Use ensemble of states over a time window to estimate temporal correlations - +# +# (work in progress) #--------------------------------------------------------------------------------------------------- def PF(self,Xb,yo): @@ -429,7 +527,7 @@ def PF(self,Xb,yo): if (Neff < edim/2): # Apply additive inflation (remove sample mean) const=1.0 - rmat=np.rand.randn(xdim,edim) * np.matlib.repmat(np.std(Xa,axis=1),1,edim) * const; + rmat=np.random.randn(xdim,edim) * np.matlib.repmat(np.std(Xa,axis=1),1,edim) * const; Xa = Xa + rmat - np.matlib.repmat(np.mean(rmat,axis=1),1,edim); KH = [0] # dummy output diff --git a/class_lorenz63.py b/class_lorenz63.py index 8137907..ccaa34a 100755 --- a/class_lorenz63.py +++ b/class_lorenz63.py @@ -290,12 +290,11 @@ def plot_lines_and_points(self,states,points,cvec,outfile='l63-3d',plot_title='L mode='lines-and-markers', marker=dict( size=1, - color=cvec, - colorscale='Viridis', + color='rgb(0,0,0)', opacity=0.5, ), line=dict( - color=cvec, + color='rgb(0,0,0)', width=1 ) ) @@ -304,10 +303,11 @@ def plot_lines_and_points(self,states,points,cvec,outfile='l63-3d',plot_title='L x=xp, y=yp, z=zp, mode='markers', marker=dict( - size=3, + size=2, color=cvec, colorscale='Viridis', - symbol='circle-open', + opacity=0.5, +# symbol='circle-open', ), # line=dict( # color=cvec, @@ -318,7 +318,7 @@ def plot_lines_and_points(self,states,points,cvec,outfile='l63-3d',plot_title='L data = [trace0, trace1] layout = dict( - width=800, + width=1200, height=700, autosize=False, title=plot_title, diff --git a/generate_analysis_3dEns.py b/generate_analysis_3dEns.py index cea9f5e..01bf37b 100644 --- a/generate_analysis_3dEns.py +++ b/generate_analysis_3dEns.py @@ -73,24 +73,42 @@ # EnKF #method='ETKF' +#----------- +# Hybrid methods +#----------- +#method='Hybrid' + #das.setMethod(method) #----------------------------------------------------------------------- # Initialize the ensemble #----------------------------------------------------------------------- xa = sv.x0 +edim = 20 #3 #20 bias_init = 0 sigma_init = 0.1 -edim = 3 #20 Xa = das.initEns(xa,mu=bias_init,sigma=sigma_init,edim=edim) +print('ensemble dimension = ') +print(edim) +print('initial bias = ') +print(bias_init) +print('initial standard deviation = ') +print(sigma_init) +print('X0 = ') +print(Xa) + +#exit() + #----------------------------------------------------------------------- # Conduct data assimilation process #----------------------------------------------------------------------- # xa = sv.x0 xa_history = np.zeros_like(x_nature) +xa_history[:] = np.nan KH_history = [] +ii=0 for i in range(0,maxit-acyc_step,acyc_step): #---------------------------------------------- @@ -128,6 +146,14 @@ xa = np.mean(Xa,axis=1).T # print('xa = ') # print(xa) + +# print('x_nature[i+acyc_step,:] = ') +# print(x_nature[i+acyc_step,:,]) + + ii=ii+1 +# if (ii>4): +# exit() + # print('KH = ') # print(KH) diff --git a/generate_analysis_4dDet.py b/generate_analysis_4dDet.py index 130e39b..39a3476 100644 --- a/generate_analysis_4dDet.py +++ b/generate_analysis_4dDet.py @@ -1 +1,2 @@ # You're on your own. +print('Good luck.') diff --git a/generate_analysis_4dEns.py b/generate_analysis_4dEns.py index 130e39b..92a65b2 100644 --- a/generate_analysis_4dEns.py +++ b/generate_analysis_4dEns.py @@ -1 +1,2 @@ # You're on your own. +print('Strength and honor') diff --git a/generate_analysis_init.py b/generate_analysis_init.py new file mode 100644 index 0000000..5c66a7e --- /dev/null +++ b/generate_analysis_init.py @@ -0,0 +1,155 @@ +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 + +#----------------------------------------------------------------------- +# Read the L63 nature run +#----------------------------------------------------------------------- +infile = 'x_nature.pkl' +sv = state_vector() +sv = sv.load(infile) +x_nature = sv.getTrajectory() +maxit,xdim = np.shape(x_nature) + +#----------------------------------------------------------------------- +# Read the L63 observations +#----------------------------------------------------------------------- +infile = 'y_obs.pkl' +obs = obs_data() +obs = obs.load(infile) + +#----------------------------------------------------------------------- +# 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 + +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,:]) + +_,ydim = np.shape(y_obs) + +#----------------------------------------------------------------------- +# Initialize the da system +#----------------------------------------------------------------------- +das = da_system() +das.setStateVector(sv) +das.setObsData(obs) +das.xdim = xdim +das.ydim = ydim +das.x0 = x_nature[0,:] +das.t = sv.getTimes() +das.t0 = das.t[0] + +#----------------------------------------------------------------------- +# Initialize the error covariances B and R, and the linearized +# observation operator H +#----------------------------------------------------------------------- + +I = np.identity(xdim) + +# Set background error covariance +sigma_b = 0.9 +B = I * sigma_b**2 + +# Set observation error covariance +sigma_r = 1.0 +R = I * sigma_r**2 + +# 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) + +das.setB(sigma_b**2*I) +das.setR(sigma_r**2*I) +das.setH(I) + +#----------------------------------------------------------------------- +# Initialize the timesteps +#----------------------------------------------------------------------- +t_nature = sv.getTimes() +acyc_step = 10 # (how frequently to perform an analysis) +dtau = (t_nature[acyc_step] - t_nature[0]) +fcst_step = acyc_step # (may need to change for 4D DA methods) +fcst_dt = dtau / fcst_step +maxit,xdim = np.shape(x_nature) + +# Store basic timing info in das object +das.acyc_step = acyc_step +das.dtau = dtau +das.fcst_step = fcst_step +das.fcst_dt = fcst_dt +das.dt = (t_nature[1] - t_nature[0]) +das.maxit = maxit +das.xdim = xdim + +#----------------------------------------------------------------------- +# Choose DA method: +#----------------------------------------------------------------------- + +#----------- +# Test basic functionality +#----------- +method='skip' + +#----------- +# 3D methods +#----------- +# Nudging +#method='nudging' +# OI +#method='OI' +# 3D-Var +#method='3DVar' + +#----------- +# Ensemble methods +#----------- +# EnKF +method='ETKF' +# Particle filter +#method='PF' + +#----------- +# 4D methods +#----------- +# 4D-Var +#method='4DVar' +# 4DEnVar +#method='4DEnVar' +# 4DETKF +#method='4DETKF' + +#----------- +# Hybrid methods +#----------- +# Hybrid-Gain +#method='Hybrid' + +das.setMethod(method) + +#----------------------------------------------------------------------- +# Store DA object +#----------------------------------------------------------------------- +name = 'x_analysis_init' +outfile=name+'.pkl' +das.save(outfile) + +print(das) diff --git a/generate_nature_run.py b/generate_nature_run.py index b497cdc..47c431a 100755 --- a/generate_nature_run.py +++ b/generate_nature_run.py @@ -14,7 +14,7 @@ t0=0.0 tf=40.0 -dt=0.01 +dt=0.001 state0 = [1.0, 1.0, 1.0] tvec = np.arange(t0, tf, dt) diff --git a/plot_analysis_vs_nature.py b/plot_analysis_vs_nature.py index 5479e91..8cdf2a7 100644 --- a/plot_analysis_vs_nature.py +++ b/plot_analysis_vs_nature.py @@ -8,7 +8,7 @@ # Read state vector objects #------------------------------------------------------------------ name = 'x_analysis' -method = 'Hybrid' +method = 'ETKF' infile1 = name+'_'+method+'.pkl' das = da_system() das = das.load(infile1) diff --git a/runall_tutorial_1.sh b/runall_tutorial_1.sh new file mode 100644 index 0000000..3864355 --- /dev/null +++ b/runall_tutorial_1.sh @@ -0,0 +1,5 @@ +python generate_nature_run.py +python generate_nature_Mhist.py +python generate_nature_QR.py +python generate_nature_LEs.py +python generate_observations.py