Skip to content

Update AnDA_analog_forecasting.py #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pnavaro wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Update AnDA_analog_forecasting.py #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
d4e2343
Update AnDA_analog_forecasting.py
pnavaro Dec 23, 2021
f564a9a
Update test_AnDA.ipynb
pnavaro Dec 23, 2021
c934f84
Update AnDA_model_forecasting.py
pnavaro Dec 23, 2021
b0f80e0
Use RandomState
pnavaro Dec 23, 2021
e413d73
put a seed in analog forecasting function was not a good idea
pnavaro Dec 25, 2021
68ec640
Update AnDA_generate_data.py
pnavaro Jan 6, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 8 additions & 8 deletions AnDA_codes/AnDA_analog_forecasting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ def AnDA_analog_forecasting(x, AF):
xf_mean = np.zeros([N,n])
stop_condition = 0
i_var = np.array([0])

# local or global analog forecasting
while (stop_condition !=1):

Expand Down Expand Up @@ -62,20 +62,20 @@ def AnDA_analog_forecasting(x, AF):
xf_tmp[:,i_var] = AF.catalog.successors[np.ix_(index_knn[i_N,:],i_var)]

# weighted mean and covariance
xf_mean[i_N,i_var] = np.sum(xf_tmp[:,i_var]*np.repeat(weights[i_N,:][np.newaxis].T,len(i_var),1),0)
E_xf = (xf_tmp[:,i_var]-np.repeat(xf_mean[i_N,i_var][np.newaxis],AF.k,0)).T
cov_xf = 1.0/(1.0-np.sum(np.power(weights[i_N,:],2)))*np.dot(np.repeat(weights[i_N,:][np.newaxis],len(i_var),0)*E_xf,E_xf.T)
xf_mean[i_N,i_var] = np.sum(xf_tmp[:,i_var]*weights[i_N,:,np.newaxis],0)
E_xf = (xf_tmp[:,i_var]-xf_mean[np.newaxis,i_N,i_var]).T
cov_xf = 1.0/(1.0-np.sum(np.power(weights[i_N,:],2)))*np.dot(weights[np.newaxis,i_N,:]*E_xf,E_xf.T)

# method "locally-incremental"
elif (AF.regression == 'increment'):

# compute the analog forecasts
xf_tmp[:,i_var] = np.repeat(x[i_N,i_var][np.newaxis],AF.k,0) + AF.catalog.successors[np.ix_(index_knn[i_N,:],i_var)]-AF.catalog.analogs[np.ix_(index_knn[i_N,:],i_var)]
xf_tmp[:,i_var] = x[np.newaxis,i_N,i_var] + AF.catalog.successors[np.ix_(index_knn[i_N,:],i_var)]-AF.catalog.analogs[np.ix_(index_knn[i_N,:],i_var)]

# weighted mean and covariance
xf_mean[i_N,i_var] = np.sum(xf_tmp[:,i_var]*np.repeat(weights[i_N,:][np.newaxis].T,len(i_var),1),0)
E_xf = (xf_tmp[:,i_var]-np.repeat(xf_mean[i_N,i_var][np.newaxis],AF.k,0)).T
cov_xf = 1.0/(1-np.sum(np.power(weights[i_N,:],2)))*np.dot(np.repeat(weights[i_N,:][np.newaxis],len(i_var),0)*E_xf,E_xf.T)
xf_mean[i_N,i_var] = np.sum(xf_tmp[:,i_var]*weights[i_N,:,np.newaxis],0)
E_xf = (xf_tmp[:,i_var]-xf_mean[np.newaxis,i_N,i_var]).T
cov_xf = 1.0/(1-np.sum(np.power(weights[i_N,:],2)))*np.dot(weights[np.newaxis,i_N,:]*E_xf,E_xf.T)

# method "locally-linear"
elif (AF.regression == 'local_linear'):
Expand Down
18 changes: 10 additions & 8 deletions AnDA_codes/AnDA_data_assimilation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,14 @@
from AnDA_codes.AnDA_stat_functions import resampleMultinomial, inv_using_SVD
from tqdm import tqdm

def AnDA_data_assimilation(yo, DA):
def AnDA_data_assimilation(yo, DA, seed=1):
"""
Apply stochastic and sequential data assimilation technics using
model forecasting or analog forecasting.
"""

rng = np.random.RandomState(seed)

# dimensions
n = len(DA.xb)
T = yo.values.shape[0]
Expand All @@ -45,7 +47,7 @@ class x_hat:
for k in tqdm(range(0,T)):
# update step (compute forecasts)
if k==0:
xf = np.random.multivariate_normal(DA.xb, DA.B, DA.N)
xf = rng.multivariate_normal(DA.xb, DA.B, DA.N)
else:
xf, m_xa_part_tmp = DA.m(x_hat.part[k-1,:,:])
m_xa_part[k,:,:] = m_xa_part_tmp
Expand All @@ -55,7 +57,7 @@ class x_hat:
# analysis step (correct forecasts with observations)
i_var_obs = np.where(~np.isnan(yo.values[k,:]))[0]
if (len(i_var_obs)>0):
eps = np.random.multivariate_normal(np.zeros(len(i_var_obs)),DA.R[np.ix_(i_var_obs,i_var_obs)],DA.N)
eps = rng.multivariate_normal(np.zeros(len(i_var_obs)),DA.R[np.ix_(i_var_obs,i_var_obs)],DA.N)
yf = np.dot(DA.H[i_var_obs,:],xf.T).T
SIGMA = np.dot(np.dot(DA.H[i_var_obs,:],Pf[k,:,:]),DA.H[i_var_obs,:].T)+DA.R[np.ix_(i_var_obs,i_var_obs)]
SIGMA_INV = np.linalg.inv(SIGMA)
Expand Down Expand Up @@ -94,7 +96,7 @@ class x_hat:
k_count = 0
m_xa_traj = []
weights_tmp = np.zeros(DA.N)
xf = np.random.multivariate_normal(DA.xb, DA.B, DA.N)
xf = rng.multivariate_normal(DA.xb, DA.B, DA.N)
i_var_obs = np.where(~np.isnan(yo.values[k,:]))[0]
if (len(i_var_obs)>0):
# weights
Expand All @@ -103,7 +105,7 @@ class x_hat:
# normalization
weights_tmp = weights_tmp/np.sum(weights_tmp)
# resampling
indic = resampleMultinomial(weights_tmp)
indic = resampleMultinomial(rng, weights_tmp)
x_hat.part[k,:,:] = xf[indic,:]
weights_tmp_indic = weights_tmp[indic]/sum(weights_tmp[indic])
x_hat.values[k,:] = sum(xf[indic,:]*weights_tmp_indic[np.newaxis].T,0)
Expand All @@ -113,12 +115,12 @@ class x_hat:
# weights
weights_tmp = 1.0/N
# resampling
indic = resampleMultinomial(weights_tmp)
indic = resampleMultinomial(rng, weights_tmp)
x_hat.weights[k,:] = weights_tmp_indic

for k in tqdm(range(1,T)):
# update step (compute forecasts) and add small Gaussian noise
xf, tej = DA.m(x_hat.part[k-1,:,:]) +np.random.multivariate_normal(np.zeros(xf.shape[1]),DA.B/100.0,xf.shape[0])
xf, tej = DA.m(x_hat.part[k-1,:,:]) + rng.multivariate_normal(np.zeros(xf.shape[1]),DA.B/100.0,xf.shape[0])
if (k_count<len(m_xa_traj)):
m_xa_traj[k_count] = xf
else:
Expand All @@ -133,7 +135,7 @@ class x_hat:
# normalization
weights_tmp = weights_tmp/np.sum(weights_tmp)
# resampling
indic = resampleMultinomial(weights_tmp)
indic = resampleMultinomial(rng, weights_tmp)
# stock results
x_hat.part[k-k_count_end:k+1,:,:] = np.asarray(m_xa_traj)[:,indic,:]
weights_tmp_indic = weights_tmp[indic]/np.sum(weights_tmp[indic])
Expand Down
23 changes: 10 additions & 13 deletions AnDA_codes/AnDA_generate_data.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from scipy.integrate import odeint
from AnDA_codes.AnDA_dynamical_models import AnDA_Lorenz_63, AnDA_Lorenz_96

def AnDA_generate_data(GD):
def AnDA_generate_data(GD, seed=1):
""" Generate the true state, noisy observations and catalog of numerical simulations. """

# initialization
Expand All @@ -29,16 +29,16 @@ class catalog:

# test on parameters
if GD.dt_states>GD.dt_obs:
print('Error: GD.dt_obs must be bigger than GD.dt_states');
raise ValueError('GD.dt_obs must be bigger than GD.dt_states');
if (np.mod(GD.dt_obs,GD.dt_states)!=0):
print('Error: GD.dt_obs must be a multiple of GD.dt_states');
raise ValueError('GD.dt_obs must be a multiple of GD.dt_states');

# use this to generate the same data for different simulations
np.random.seed(1);
rng = np.random.RandomState(seed);

if (GD.model == 'Lorenz_63'):

# 5 time steps (to be in the attractor space)
# Go to time = 5 to be in the attractor space
x0 = np.array([8.0,0.0,30.0]);
S = odeint(AnDA_Lorenz_63,x0,np.arange(0,5+0.000001,GD.dt_integration),args=(GD.parameters.sigma,GD.parameters.rho,GD.parameters.beta));
x0 = S[S.shape[0]-1,:];
Expand All @@ -51,7 +51,7 @@ class catalog:
xt.values = S[t_xt,:];

# generate partial/noisy observations (yo)
eps = np.random.multivariate_normal(np.zeros(3),GD.sigma2_obs*np.eye(3,3),T_test);
eps = rng.multivariate_normal(np.zeros(3),GD.sigma2_obs*np.eye(3,3),T_test);
yo_tmp = S[t_xt,:]+eps[t_xt,:];
t_yo = np.arange(0,T_test,GD.dt_obs);
i_t_obs = np.where((np.in1d(t_xt,t_yo))==True)[0];
Expand All @@ -62,15 +62,15 @@ class catalog:
#generate catalog
S = odeint(AnDA_Lorenz_63,S[S.shape[0]-1,:],np.arange(0.01,GD.nb_loop_train+0.000001,GD.dt_integration),args=(GD.parameters.sigma,GD.parameters.rho,GD.parameters.beta));
T_train = S.shape[0];
eta = np.random.multivariate_normal(np.zeros(3),GD.sigma2_catalog*np.eye(3,3),T_train);
eta = rng.multivariate_normal(np.zeros(3),GD.sigma2_catalog*np.eye(3,3),T_train);
catalog_tmp = S+eta;
catalog.analogs = catalog_tmp[0:-GD.dt_states,:];
catalog.successors = catalog_tmp[GD.dt_states:,:]
catalog.source = GD.parameters;

elif (GD.model == 'Lorenz_96'):

# 5 time steps (to be in the attractor space)
# Go to time = 5 to be in the attractor space
x0 = GD.parameters.F*np.ones(GD.parameters.J);
x0[np.int(np.around(GD.parameters.J/2))] = x0[np.int(np.around(GD.parameters.J/2))] + 0.01;
S = odeint(AnDA_Lorenz_96,x0,np.arange(0,5+0.000001,GD.dt_integration),args=(GD.parameters.F,GD.parameters.J));
Expand All @@ -84,7 +84,7 @@ class catalog:
xt.values = S[t_xt,:];

# generate partial/noisy observations (yo)
eps = np.random.multivariate_normal(np.zeros(GD.parameters.J),GD.sigma2_obs*np.eye(GD.parameters.J),T_test);
eps = rng.multivariate_normal(np.zeros(GD.parameters.J),GD.sigma2_obs*np.eye(GD.parameters.J),T_test);
yo_tmp = S[t_xt,:]+eps[t_xt,:];
t_yo = np.arange(0,T_test,GD.dt_obs);
i_t_obs = np.where((np.in1d(t_xt,t_yo))==True)[0];
Expand All @@ -95,15 +95,12 @@ class catalog:
# generate catalog
S = odeint(AnDA_Lorenz_96,S[S.shape[0]-1,:],np.arange(0.01,GD.nb_loop_train+0.000001,GD.dt_integration),args=(GD.parameters.F,GD.parameters.J));
T_train = S.shape[0];
eta = np.random.multivariate_normal(np.zeros(GD.parameters.J),GD.sigma2_catalog*np.eye(GD.parameters.J,GD.parameters.J),T_train);
eta = rng.multivariate_normal(np.zeros(GD.parameters.J),GD.sigma2_catalog*np.eye(GD.parameters.J,GD.parameters.J),T_train);
catalog_tmp = S+eta;
catalog.analogs = catalog_tmp[0:-GD.dt_states,:];
catalog.successors = catalog_tmp[GD.dt_states:,:]
catalog.source = GD.parameters;

# reinitialize random generator number
np.random.seed()

return catalog, xt, yo;


Expand Down
19 changes: 9 additions & 10 deletions AnDA_codes/AnDA_model_forecasting.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,19 +18,18 @@ def AnDA_model_forecasting(x,GD):
# initializations
N, n = x.shape;
xf = np.zeros([N,n]);
xf_mean = np.zeros([N,n]);

if (GD.model == 'Lorenz_63'):
for i_N in range(0,N):
S = odeint(AnDA_Lorenz_63,x[i_N,:],np.arange(0,GD.dt_integration+0.000001,GD.dt_integration),args=(GD.parameters.sigma,GD.parameters.rho,GD.parameters.beta));
xf[i_N,:] = S[-1,:];
fmodel = {}

elif (GD.model == 'Lorenz_96'):
for i_N in range(0,N):
S = odeint(AnDA_Lorenz_96,x[i_N,:],np.arange(0,GD.dt_integration+0.000001,GD.dt_integration),args=(GD.parameters.F,GD.parameters.J));
xf[i_N,:] = S[-1,:];
fmodel['Lorenz_63'] = lambda x0: odeint(AnDA_Lorenz_63,x0,(0,GD.dt_integration),
args=(GD.parameters.sigma,GD.parameters.rho,GD.parameters.beta))[-1]

fmodel['Lorenz_96'] = lambda x0: odeint(AnDA_Lorenz_96,x0,(0,GD.dt_integration),
args=(GD.parameters.F,GD.parameters.J))[-1]

xf = np.apply_along_axis(fmodel[GD.model], 1, x)
xf_mean = np.copy(xf)

xf_mean = xf;
return xf, xf_mean


4 changes: 2 additions & 2 deletions AnDA_codes/AnDA_stat_functions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ def sample_discrete(prob, r, c):
M = M+1*(R>cumprob[i]);
return int(M)

def resampleMultinomial(w):
def resampleMultinomial(rng, w):
""" Multinomial resampler. """

M = np.max(w.shape);
Expand All @@ -62,7 +62,7 @@ def resampleMultinomial(w):
i = 0;
indx = [];
while (i<=(M-1)):
sampl = np.random.rand(1,1);
sampl = rng.rand(1,1);
j = 0;
while (Q[j]<sampl):
j = j+1;
Expand Down
3 changes: 1 addition & 2 deletions test_AnDA.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -286,8 +286,7 @@
"### COMPARISON BETWEEN CLASSICAL AND ANALOG DATA ASSIMILATION\n",
"\n",
"# plot\n",
"fig=plt.figure()\n",
"ax=fig.gca(projection='3d')\n",
"ax = plt.figure().add_subplot(projection='3d')\n",
"line1,=ax.plot(xt.values[:,0],xt.values[:,1],xt.values[:,2],'k')\n",
"line2,=ax.plot(x_hat_analog.values[:,0],x_hat_analog.values[:,1],x_hat_analog.values[:,2],'r')\n",
"line3,=ax.plot(x_hat_classical.values[:,0],x_hat_classical.values[:,1],x_hat_classical.values[:,2],'b')\n",
Expand Down