<jupyter_start><jupyter_text>Dataset generation<jupyter_code>## second u0 is off-phase from the first by half a Lyapunov time
u0 = np.array([7.432487609628195, 10.02071718705213, 29.62297428638419])

dt        = 0.01
t_lyap    = 0.9**(-1)
N_lyap    = int(t_lyap/dt)  # number of time steps in one Lyapunov time

# number of time steps for washout, train, validation and PH window
N_washout = 1*N_lyap
N_train   = 8*N_lyap
N_val     = 3*N_lyap
N_tstart  = 24*N_lyap #start for the test set
N_test    = 500*N_lyap

norm      = np.array([35.23020746, 45.09776766, 36.07598481])

# generate data 
U         = solve_ode(N_washout+N_train+N_val+N_test, dt, u0)

# washout
U_washout = U[:N_washout]
# training
U_t       = U[N_washout:N_washout+N_train-1]
Y_t       = U[N_washout+1:N_washout+N_train]
# training + validation
U_tv      = U[N_washout:N_washout+N_train+N_val-1]
Y_tv      = U[N_washout+1:N_washout+N_train+N_val]
# validation
Y_val     = U[N_washout+N_train:N_washout+N_train+N_val]<jupyter_output><empty_output><jupyter_text>Import Results From  Optimization Runs<jupyter_code>model_informed = True
if model_informed:
    string = '_MI.h5'
    %run ./Functions_MI.ipynb
else:
    string = '.h5'

#BO and Grid Search in SSV
hf       = h5py.File('./data/Lor_short_SSV_50_5' + str(string),'r')
Min_25G  = np.array(hf.get('minimum'))
hf.close()

hf       = h5py.File('./data/Lor_short_SSV_50_7' + str(string),'r')
Min_50G  = np.array(hf.get('minimum'))
hf.close()

# #BO and Grid Search in RVC 
hf       = h5py.File('./data/Lor_short_RVC_50_5' + str(string),'r')
Min_25G_mmv  = np.array(hf.get('minimum'))
hf.close()

hf       = h5py.File('./data/Lor_short_RVC_50_7' + str(string),'r')
Min_50G_mmv  = np.array(hf.get('minimum'))
hf.close()

# BO and Grid Search in WFV
hf       = h5py.File('./data/Lor_short_WFV_50_5' + str(string),'r')
Min_25G_wfv  = np.array(hf.get('minimum'))
hf.close()

hf       = h5py.File('./data/Lor_short_WFV_50_7' + str(string),'r')
Min_50G_wfv  = np.array(hf.get('minimum'))
hf.close()

#BO and Grid Search in KFV 
hf       = h5py.File('./data/Lor_short_KFV_50_5' + str(string),'r')
Min_25G_kfv  = np.array(hf.get('minimum'))
hf.close()

hf       = h5py.File('./data/Lor_short_KFV_50_7' + str(string),'r')
Min_50G_kfv  = np.array(hf.get('minimum'))
hf.close()

#BO and Grid Search in RV 
hf       = h5py.File('./data/Lor_short_RV_50_5' + str(string),'r')
Min_25G_mv  = np.array(hf.get('minimum'))
hf.close()

hf       = h5py.File('./data/Lor_short_RV_50_7' + str(string),'r')
Min_50G_mv  = np.array(hf.get('minimum'))
hf.close()

#BO and Grid Search in KFC
hf       = h5py.File('./data/Lor_short_KFC_50_5' + str(string),'r')
Min_25G_kfo  = np.array(hf.get('minimum'))
hf.close()

hf       = h5py.File('./data/Lor_short_KFC_50_7' + str(string),'r')
Min_50G_kfo  = np.array(hf.get('minimum'))
hf.close()

#BO and Grid Search WFC
hf       = h5py.File('./data/Lor_short_WFC_50_5' + str(string),'r')
Min_25G_wfc  = np.array(hf.get('minimum'))
hf.close()

hf       = h5py.File('./data/Lor_short_WFC_50_7' + str(string),'r')
Min_50G_wfc  = np.array(hf.get('minimum'))
hf.close()<jupyter_output><empty_output><jupyter_text>ESN Initiliazation Parameters<jupyter_code>bias_in = 1.0 #input bias
bias_out = 1.0 #output bias 
N_units = 100 #units in the reservoir
dim = 3 # dimension of inputs (and outputs) 
connectivity   = 3 

sigma_in = 1.0 #input scaling
rho = 1.0  # spectral radius
sparseness =  1 - connectivity/(N_units-1) 

tikh = 1e-11  # Tikhonov factor
# Functions in the test set

def f(x):
    """Computes MSE and PH in the test set of N_test points"""
    
    global rho, sigma_in
    rho      = x[0]
    sigma_in = x[1] 
                
    #Train on the entire dataset
    Wout     = train(U_washout, U_tv, Y_tv, tikh)[1]

    N_test   = 100 # N_test intervals
    Mean_MSE = 0
    Mean_PH  = 0
    kk       = 5   # the PH is computed up to kk*N_val for each interval
    
    #Different Folds in the cross-validation
    for i in range(N_test):
        
        # data for washout and target in each interval
        U_wash    = U[N_tstart - N_washout +i*N_val : N_tstart+i*N_val]
        Y_t_MSE   = U[N_tstart+i*N_val:N_tstart+i*N_val+N_val] 
        Y_t_PH    = U[N_tstart+i*N_val:N_tstart+i*N_val+kk*N_val] 
        
        #washout for each interval
        xa1       = open_loop(U_wash, np.zeros(N_units))[-1]
        
        # Mean Square Error
        Yh_t_MSE  = closed_loop(N_val-1, xa1, Wout)[0]
        Mean1     = np.log10(np.mean(((Yh_t_MSE-Y_t_MSE))**2))
        # np.isnan because trajectory may become unbounded in model-informed architecture
        if np.isnan(Mean1):
            Mean1 = 10
        Mean_MSE += Mean1
        
        # Prediction Horizon
        Mean_PH  += predictability_horizon(xa1,Y_t_PH,Wout)
                
    return Mean_MSE/N_test, Mean_PH/N_test<jupyter_output><empty_output><jupyter_text>Run Ensemble in the test set<jupyter_code>%%time

#Compute Ensemble
ensemble = 50

# Initialize the arrays 
res_Gr     = np.zeros((ensemble,2))
res_BO     = np.zeros((ensemble,2))
res_Gr_mmv = np.zeros((ensemble,2))
res_BO_mmv = np.zeros((ensemble,2))
res_Gr_wfv = np.zeros((ensemble,2))
res_BO_wfv = np.zeros((ensemble,2))
res_Gr_kfv = np.zeros((ensemble,2))
res_BO_kfv = np.zeros((ensemble,2))
res_Gr_mv   = np.zeros((ensemble,2))
res_BO_mv   = np.zeros((ensemble,2))
res_Gr_kfo  = np.zeros((ensemble,2))
res_BO_kfo  = np.zeros((ensemble,2))
res_Gr_wfc  = np.zeros((ensemble,2))
res_BO_wfc  = np.zeros((ensemble,2))

for i in range(ensemble):
     
    print('Ensemble    :',i+1)
        
    # Win and W generation
    seed= i+1
    rnd = np.random.RandomState(seed)

    Win = np.zeros((dim+1, N_units))
    for j in range(N_units):
        Win[rnd.randint(0, dim+1),j] = rnd.uniform(-1, 1) #only one element different from zero per row
    
    # practical way to set the sparseness
    W = rnd.uniform(-1, 1, (N_units, N_units)) * (rnd.rand(N_units, N_units) < (1-sparseness))
    spectral_radius = np.max(np.abs(np.linalg.eigvals(W)))
    W /= spectral_radius #scaled to have unitary spec radius
    
    #Compute the performance in the test set for the best hyperparameters
    temp = f(Min_50G[i,:2])
    res_Gr[i,0]  = temp[0]
    res_Gr[i,1]  = temp[1]
        
    temp = f(Min_25G[i,:2])
    res_BO[i,0]  = temp[0]
    res_BO[i,1]  = temp[1]
        
    temp = f(Min_50G_mmv[i,:2])
    res_Gr_mmv[i,0]  = temp[0]
    res_Gr_mmv[i,1]  = temp[1]
    
    temp = f(Min_25G_mmv[i,:2])
    res_BO_mmv[i,0]  = temp[0]
    res_BO_mmv[i,1]  = temp[1]
        
    temp = f(Min_50G_wfv[i,:2])
    res_Gr_wfv[i,0]  = temp[0]
    res_Gr_wfv[i,1]  = temp[1]
    
    temp = f(Min_25G_wfv[i,:2])
    res_BO_wfv[i,0]  = temp[0]
    res_BO_wfv[i,1]  = temp[1]
    
    temp = f(Min_50G_kfv[i,:2])
    res_Gr_kfv[i,0]  = temp[0]
    res_Gr_kfv[i,1]  = temp[1]
    
    temp = f(Min_25G_kfv[i,:2])
    res_BO_kfv[i,0]  = temp[0]
    res_BO_kfv[i,1]  = temp[1]
    
    temp = f(Min_50G_kfo[i,:2])
    res_Gr_kfo[i,0]  = temp[0]
    res_Gr_kfo[i,1]  = temp[1]
    
    temp = f(Min_25G_kfo[i,:2])
    res_BO_kfo[i,0]  = temp[0]
    res_BO_kfo[i,1]  = temp[1]
            
    temp = f(Min_50G_mv[i,:2])
    res_Gr_mv[i,0]  = temp[0]
    res_Gr_mv[i,1]  = temp[1]
    
    temp = f(Min_25G_mv[i,:2])
    res_BO_mv[i,0]  = temp[0]
    res_BO_mv[i,1]  = temp[1]

    temp = f(Min_50G_wfc[i,:2])
    res_Gr_wfc[i,0]  = temp[0]
    res_Gr_wfc[i,1]  = temp[1]
    
    temp = f(Min_25G_wfc[i,:2])
    res_BO_wfc[i,0]  = temp[0]
    res_BO_wfc[i,1]  = temp[1]

#Save Results

fln = './data/Short_Lorenz_PostProc' + str(string)

hf = h5py.File(fln,'w')
hf.create_dataset('res_Gr     ',data=res_Gr)
hf.create_dataset('res_BO     ',data=res_BO)
hf.create_dataset('res_Gr_mmv ',data=res_Gr_mmv)
hf.create_dataset('res_BO_mmv ',data=res_BO_mmv)
hf.create_dataset('res_Gr_mv  ',data=res_Gr_mv)
hf.create_dataset('res_BO_mv  ',data=res_BO_mv)
hf.create_dataset('res_Gr_kfv ',data=res_Gr_kfv)
hf.create_dataset('res_BO_kfv ',data=res_BO_kfv)
hf.create_dataset('res_Gr_kfo ',data=res_Gr_kfo)
hf.create_dataset('res_BO_kfo ',data=res_BO_kfo)
hf.create_dataset('res_Gr_wfv ',data=res_Gr_wfv)
hf.create_dataset('res_BO_wfv ',data=res_BO_wfv)
hf.create_dataset('res_Gr_wfc ',data=res_Gr_wfc)
hf.create_dataset('res_BO_wfc ',data=res_BO_wfc)
hf.close()<jupyter_output><empty_output><jupyter_text>Run Ensmeble in the test set with fixed hyperparametersWe compute the performance of the ensemble using the same hyperparameters for all the networks.The hyperparameters used are the ones obtained through Bayesian Optimization in the first network using chaotic Recycle Validation and chaotic K-Fold Validation.<jupyter_code>%%time

model_informed = False
if model_informed:
    %run ./Functions_MI.ipynb

#Compute Ensemble
ensemble = 50

# Initialize the arrays 
res_BO_mmv_fix = np.zeros((ensemble,2))
res_BO_kfo_fix = np.zeros((ensemble,2))

for i in range(ensemble):
     
    print('Ensemble    :',i+1)
        
    # Win and W generation
    seed= i+1
    rnd = np.random.RandomState(seed)

    Win = np.zeros((dim+1, N_units))
    for j in range(N_units):
        Win[rnd.randint(0, dim+1),j] = rnd.uniform(-1, 1) #only one element different from zero per row
    
    # practical way to set the sparseness
    W = rnd.uniform(-1, 1, (N_units, N_units)) * (rnd.rand(N_units, N_units) < (1-sparseness))
    spectral_radius = np.max(np.abs(np.linalg.eigvals(W)))
    W /= spectral_radius #scaled to have unitary spec radius
    
    #Compute the performance in the test set for the best hyperparameters
    temp = f(Min_25G_mmv[0,:2])
    res_BO_mmv_fix[i,0]  = temp[0]
    res_BO_mmv_fix[i,1]  = temp[1]
        
    temp = f(Min_25G_kfo[0,:2])
    res_BO_kfo_fix[i,0]  = temp[0]
    res_BO_kfo_fix[i,1]  = temp[1]

#Save Results

fln = './data/Short_Lorenz_PostProc_fix.h5'

hf = h5py.File(fln,'w')
hf.create_dataset('res_BO_mmv_fix',data=res_BO_mmv_fix)
hf.create_dataset('res_BO_kfo_fix',data=res_BO_kfo_fix)
hf.close()<jupyter_output><empty_output><jupyter_text>Comparison SSV in validation and 12 to 15 LTs<jupyter_code>%%time

from skopt.plots import plot_convergence

# We use only one interval for the test set (12 to 15LTs) to compute the MSE
# Test1 is implemented in Functions.ipynb

#Compute Ensemble
ensemble = 50

#Computing 
res_Gr1     = np.zeros((ensemble,1))
res_BO1     = np.zeros((ensemble,1))

for i in range(ensemble):
        
    # Win and W generation
    seed =i+1
    rnd = np.random.RandomState(seed)

    Win = np.zeros((dim+1, N_units))
    for j in range(N_units):
        Win[rnd.randint(0, dim+1),j] = rnd.uniform(-1, 1)
        
    W = rnd.uniform(-1, 1, (N_units, N_units)) * (rnd.rand(N_units, N_units) < (1-sparseness))
    spectral_radius = np.max(np.abs(np.linalg.eigvals(W)))
    W /= spectral_radius
    
    # 
    temp = Test1(Min_50G[i,:2])
    res_Gr1[i,0]  = 10**temp
    
    temp = Test1(Min_25G[i,:2])
    res_BO1[i,0]  = 10**temp
    
fln = './data/Short_Lorenz_PostProc2.h5'

hf = h5py.File(fln,'w')
hf.create_dataset('res_Gr1     ',data=res_Gr1)
hf.create_dataset('res_BO1     ',data=res_BO1)
hf.create_dataset('Min_25G     ',data=10**Min_25G)
hf.create_dataset('Min_50G     ',data=10**Min_50G)
hf.close()
## Model-informed<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1442905345409341, 'fraction_numerical': 0.03933581957023619, 'mean_word_length': 3.5232931726907633, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '256942', 'license': 'MIT', 'n_tokens_mistral': 5182, 'n_tokens_neox': 4919, 'n_words': 858, 'path': 'source/Test_Short_Lorenz.ipynb', 'repo_name': 'MagriLab/Robust-Validation-ESN', 'size': '11262'}