From 6e9c253323d68ed7c83924d4ad935cc3cdc93a6d Mon Sep 17 00:00:00 2001 From: Sebastian Schubert Date: Fri, 22 Sep 2017 16:28:16 +0200 Subject: [PATCH] order memmap --- lorenz96_compute_secondary_instabilities.py | 24 ++++---- lorenz96_secondary_instabilities.py | 28 ++++----- lorenz96_test_secondary_instabilities.py | 13 ++-- plot_results_secondary_instabilities.py | 24 ++++---- ...lts_secondary_instabilities_projections.py | 17 +++--- plot_test_secondary_instabilities.py | 60 +++++++++---------- 6 files changed, 82 insertions(+), 84 deletions(-) diff --git a/lorenz96_compute_secondary_instabilities.py b/lorenz96_compute_secondary_instabilities.py index 1ad11d1..471d1ef 100644 --- a/lorenz96_compute_secondary_instabilities.py +++ b/lorenz96_compute_secondary_instabilities.py @@ -32,14 +32,14 @@ 'dimY' : 10, 'RescaledY' : False, 'expname' : 'secondaryinstabilities_1layer', - 'time' : np.arange(0,500,0.01), + 'time' : np.arange(0,500,0.1), 'spinup' : 100, '2lay' : False } testzeroclv=True -steplengthforsecondorder = np.arange(0,100,1) +steplengthforsecondorder = np.arange(0,200,1) hs=[1.0]#, 0.5] # , 0.0625, 0.125 , 0.25 , 0.5 , 1. ] precision='float64' for paraL96,h in product([paraL96_1lay],hs): @@ -62,18 +62,18 @@ t = paraL96['time'] - CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64') - lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(t),M),dtype='float64') - invCLV = np.memmap(savename+'/invCLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64') + CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64', order = 'F') + lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(t),M),dtype='float64', order = 'F') + invCLV = np.memmap(savename+'/invCLV.dat',mode='r',shape=(len(t),dimN,M),dtype='float64', order = 'F') lensteparray = len(steplengthforsecondorder) - contracted_CLVs = np.memmap(savename+'/contracted_clvs.dat',mode='w+',shape=(len(t),dimN,M),dtype=precision) #(final_clv,init_clv) - solution = np.memmap(savename+'/solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision) - full_solution = np.memmap(savename+'/full_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision) - #normalized_solution = np.memmap(savename+'/normalized_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision) - growth = np.memmap(savename+'/growth.dat',mode='w+',shape=(len(t),lensteparray,M),dtype=precision) + contracted_CLVs = np.memmap(savename+'/contracted_clvs.dat',mode='w+',shape=(len(t),dimN,M),dtype=precision, order = 'F') #(final_clv,init_clv) + solution = np.memmap(savename+'/solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision, order = 'F') + full_solution = np.memmap(savename+'/full_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision, order = 'F') + #normalized_solution = np.memmap(savename+'/normalized_solution.dat',mode='w+',shape=(len(t),lensteparray,dimN,M),dtype=precision, order = 'F') + growth = np.memmap(savename+'/growth.dat',mode='w+',shape=(len(t),lensteparray,M),dtype=precision, order = 'F') for tn, (ts,te) in enumerate(zip(t[0:-2],t[1:-1])): print(tn) contracted_CLVs[tn,:,:]=1/2*contract_func(invCLV[tn,:,:],CLV[tn,:,:]) @@ -92,13 +92,13 @@ dtau * np.multiply(np.exp(-dtau*np.memmap.sum(lyaploc_clv[tni:tn+n_step,:,np.newaxis], axis =0)),np.multiply(contracted_CLVs[tn+n_step,:,:],np.exp(2.0*np.memmap.sum(dtau*lyaploc_clv[tni:tn+n_step,np.newaxis,:], axis =0)))) ) for n_ind, n_step in enumerate(steplengthforsecondorder): - solution[tn,n_ind,:,:] = np.multiply(solution[tn,n_ind,:,:] ,np.exp(np.sum(dtau*lyaploc_clv[tn:tn+n_step,:,np.newaxis], axis =0))) + solution[tn,n_ind,:,:] = np.multiply(solution[tn,n_ind,:,:] ,np.exp(np.sum(dtau*lyaploc_clv[tn:tn+n_step,:,np.newaxis], axis = 0))) if n_ind == 0: continue else: full_solution[tn,n_ind,:,:] = np.matmul(CLV[tn+n_step,:,:],solution[tn,n_step,:,:]) growth[tn,n_ind,:]=np.log(np.linalg.norm(full_solution[tn,n_ind,:,:],axis=0))/(n_step*dtau) - if tn % 400 == 0: + if tn % 1 == 0: np.memmap.flush(solution) np.memmap.flush(contracted_CLVs) #np.memmap.flush(normalized_solution) diff --git a/lorenz96_secondary_instabilities.py b/lorenz96_secondary_instabilities.py index b8546bb..f1d24fa 100644 --- a/lorenz96_secondary_instabilities.py +++ b/lorenz96_secondary_instabilities.py @@ -34,7 +34,7 @@ 'dimY' : 10, 'RescaledY' : False, 'expname' : 'secondaryinstabilities_1layer', - 'time' : np.arange(0,500,0.01), + 'time' : np.arange(0,500,0.1), 'spinup' : 100, '2lay' : False } @@ -56,7 +56,7 @@ dimN = paraL96['dimX'] t = paraL96['time'] - dt = np.mean(np.diff(t)) + dt = 0.01#np.mean(np.diff(t)) savename=paraL96['expname']+"_h_"+str(h) spinup = paraL96['spinup'] @@ -66,17 +66,17 @@ if not os.path.exists(savename): os.mkdir(savename) - CLV = np.memmap(savename+'\CLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64') - BLV = np.memmap(savename+'\BLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64') - R = np.memmap(savename+'\R.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64') - lyapmean_blv = np.memmap(savename+'\lyapmean_blv.dat',mode='w+',shape=(M),dtype='float64') - lyapmean_clv = np.memmap(savename+'\lyapmean_clv.dat',mode='w+',shape=(M),dtype='float64') - lyaploc_clv = np.memmap(savename+'\lyaploc_clv',mode='w+',shape=(len(t),M),dtype='float64') - lyaploc_blv = np.memmap(savename+'\lyaploc_blv',mode='w+',shape=(len(t)-1,M),dtype='float64') - np.save(savename+'\t',t) - trajectory = np.memmap(savename+'\trajectory.dat',mode='w+',shape=(len(t),dimN),dtype='float64') - if testzeroclv: tendency = np.memmap(savename+'\tendency.dat',mode='w+',shape=(len(t),dimN),dtype='float64') - if testzeroclv: tendcorr = np.memmap(savename+'\tendcorr.dat',mode='w+',shape=(len(t)),dtype='float64') + CLV = np.memmap(savename+'/CLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64', order = 'F') + BLV = np.memmap(savename+'/BLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64', order = 'F') + R = np.memmap(savename+'/R.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64', order = 'F') + lyapmean_blv = np.memmap(savename+'/lyapmean_blv.dat',mode='w+',shape=(M),dtype='float64', order = 'F') + lyapmean_clv = np.memmap(savename+'/lyapmean_clv.dat',mode='w+',shape=(M),dtype='float64', order = 'F') + lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='w+',shape=(len(t),M),dtype='float64', order = 'F') + lyaploc_blv = np.memmap(savename+'/lyaploc_blv',mode='w+',shape=(len(t)-1,M),dtype='float64', order = 'F') + np.save(savename+'/t',t) + trajectory = np.memmap(savename+'/trajectory.dat',mode='w+',shape=(len(t),dimN),dtype='float64', order = 'F') + if testzeroclv: tendency = np.memmap(savename+'/tendency.dat',mode='w+',shape=(len(t),dimN),dtype='float64', order = 'F') + if testzeroclv: tendcorr = np.memmap(savename+'/tendcorr.dat',mode='w+',shape=(len(t)),dtype='float64', order = 'F') # Compute the exponents @@ -164,7 +164,7 @@ np.save(savename+"/paraL96",paraL96) np.save(savename+"/h",h) - invCLV = np.memmap(savename+'\invCLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64') + invCLV = np.memmap(savename+'/invCLV.dat',mode='w+',shape=(len(t),dimN,M),dtype='float64') for tn, (ts,te) in enumerate(zip(t[0:-1],t[1:])): invCLV[tn,:,:]=np.linalg.inv(CLV[tn,:,:]) diff --git a/lorenz96_test_secondary_instabilities.py b/lorenz96_test_secondary_instabilities.py index b4de784..96b121f 100644 --- a/lorenz96_test_secondary_instabilities.py +++ b/lorenz96_test_secondary_instabilities.py @@ -55,19 +55,19 @@ integrator = 'rk4' # first test clv -epsilons=10.0**np.arange(-8,1.5,0.1,dtype = precision) -intsteps=np.arange(0,100,1) +epsilons=10.0**np.arange(-2,1.6,0.1,dtype = precision) +intsteps=np.arange(0,200,1) CLVs=[1,8,14,30] # ,13, 30]#np.arange(1,2,1) -timeintervall = range(2000,4000,500) +timeintervall = range(2000,4000+1,300) for paraL96,h in product(experiments ,hs): if not paraL96['2lay'] and not h == 1.0: print("1 lay only with h = 1.");break savename=paraL96['expname']+"_h_"+str(h) spinup = paraL96['spinup'] paraL96=np.load(savename+"/paraL96.npy") paraL96=paraL96[()] - dt = np.mean(np.diff(paraL96['time'])) + dt = 0.01#np.mean(np.diff(paraL96['time'])) # M number exponents if paraL96['2lay']: M = paraL96['dimX'] + paraL96['dimX']*paraL96['dimY'] # -1 full spectrum @@ -123,12 +123,14 @@ print(' eps: '+str(epsilon)) for nc, clv in enumerate(CLVs): field.x['back']=trajectory[tn,:]+epsilon*CLV[tn,:,clv-1] + #print(0,np.log10(epsilon),np.log10(np.linalg.norm(epsilon*CLV[tn,:,clv-1]))) + for step, (a,b) in enumerate(zip(intsteps[0:-1],intsteps[1:])): stepsize=b #print('g') for counting in np.arange(0,b-a): field.integrate_back(dtau,dt =dt) - #print('h') + #print('counting') nonlinear_prediction = field.x['back'] - trajectory[tn+stepsize,:] firstorder_prediction = CLV[tn+stepsize,:,clv-1]*np.exp(dtau*np.memmap.sum(lyaploc_clv[tn:tn+stepsize,clv-1]))*epsilon secondorder_prediction = full_solution[tn,stepsize,:,clv-1]*epsilon**2.0 @@ -141,6 +143,7 @@ +CLV[tn+stepsize,:,clv-1] *np.exp(dtau*np.memmap.sum(lyaploc_clv[tn:tn+stepsize,clv-1]))) + #print(step+1,np.log10(epsilon),np.log10(np.linalg.norm(nonlinear_prediction - firstorder_prediction)),np.log10(np.linalg.norm(nonlinear_prediction - firstandsecondorder_prediction))) correlation[nc][i,en,step] = corr(nonlinear_prediction,firstorder_prediction) diff --git a/plot_results_secondary_instabilities.py b/plot_results_secondary_instabilities.py index 55d79cd..e27d44b 100644 --- a/plot_results_secondary_instabilities.py +++ b/plot_results_secondary_instabilities.py @@ -11,20 +11,20 @@ -CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64') -BLV = np.memmap(savename+'/BLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64') -R = np.memmap(savename+'/R.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64') -lyapmean_blv = np.memmap(savename+'/lyapmean_blv.dat',mode='r',shape=(M,1),dtype='float64') -lyapmean_clv = np.memmap(savename+'/lyapmean_clv.dat',mode='r',shape=(M,1),dtype='float64') -lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(t),M,1),dtype='float64') -lyaploc_blv = np.memmap(savename+'/lyaploc_blv',mode='r',shape=(len(t)-1,M,1),dtype='float64') +CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64', order = 'F') +BLV = np.memmap(savename+'/BLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64', order = 'F') +R = np.memmap(savename+'/R.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64', order = 'F') +lyapmean_blv = np.memmap(savename+'/lyapmean_blv.dat',mode='r',shape=(M,1),dtype='float64', order = 'F') +lyapmean_clv = np.memmap(savename+'/lyapmean_clv.dat',mode='r',shape=(M,1),dtype='float64', order = 'F') +lyaploc_clv = np.memmap(savename+'/lyaploc_clv',mode='r',shape=(len(t),M,1),dtype='float64', order = 'F') +lyaploc_blv = np.memmap(savename+'/lyaploc_blv',mode='r',shape=(len(t)-1,M,1),dtype='float64', order = 'F') precision='float32' -contracted_CLVs = np.memmap(savename+'/contracted_clvs.dat',mode='r',shape=(len(t),dimN,M),dtype=precision) #(final_clv,init_clv) -solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),dimN,M),dtype=precision) -normalized_solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),dimN,M),dtype=precision) -growth = np.memmap(savename+'/growth.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),M),dtype=precision) -invCLV = np.memmap(savename+'/invCLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64') +contracted_CLVs = np.memmap(savename+'/contracted_clvs.dat',mode='r',shape=(len(t),dimN,M),dtype='float64', order = 'F') #(final_clv,init_clv) +solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),dimN,M),dtype='float64', order = 'F') +normalized_solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),dimN,M),dtype='float64', order = 'F') +growth = np.memmap(savename+'/growth.dat',mode='r',shape=(len(steplengthforsecondorder),len(t),M),dtype='float64', order = 'F') +invCLV = np.memmap(savename+'/invCLV.dat',mode='r',shape=(len(t),dimN,M,1),dtype='float64', order = 'F') diff --git a/plot_results_secondary_instabilities_projections.py b/plot_results_secondary_instabilities_projections.py index 6bb0ebb..cc3d694 100644 --- a/plot_results_secondary_instabilities_projections.py +++ b/plot_results_secondary_instabilities_projections.py @@ -31,7 +31,7 @@ 'dimY' : 10, 'RescaledY' : False, 'expname' : 'secondaryinstabilities_1layer', - 'time' : np.arange(0,1000,0.1), + 'time' : np.arange(0,100,0.1), 'spinup' : 100, '2lay' : False } @@ -52,7 +52,7 @@ compute = True #compute array projections -averageintervall=np.arange(2000,8000) +averageintervall=np.arange(1000,4000) for paraL96,h in product(experiments,hs): @@ -74,8 +74,8 @@ CLV = np.memmap(savename+'/CLV.dat',mode='r',shape=(len(paraL96['time']),dimN,M),dtype='float64') if compute: - solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype=precision) - full_solution = np.memmap(savename+'/full_solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype=precision) + solution = np.memmap(savename+'/solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype=precision,order = 'F') + full_solution = np.memmap(savename+'/full_solution.dat',mode='r',shape=(len(paraL96['time']),len(steplengthforsecondorder),dimN,M),dtype=precision, order = 'F') if compute: projections=np.zeros((len(steplengthforsecondorder),dimN,M)) correlations=np.zeros((len(steplengthforsecondorder),dimN,M)) @@ -84,10 +84,11 @@ correlations = np.load(savename+"/correlations.npy") if compute: dummy = np.zeros((len(steplengthforsecondorder),solution.shape[2],solution.shape[3])) - for tn in averageintervall: - dummy = dummy + np.abs(solution[tn,:,:,:]) - dummy = dummy/len(averageintervall) - +# for tn in averageintervall: +# print(tn) +# dummy = dummy + np.abs(solution[tn,:,:,:]) +# dummy = dummy/len(averageintervall) + dummy = np.memmap.mean(np.abs(solution[averageintervall,:,:,:])) projections[1:,:,:]=np.divide(dummy[1:,:,:], np.linalg.norm(dummy[1:,:,:], axis = 1, keepdims=True)) diff --git a/plot_test_secondary_instabilities.py b/plot_test_secondary_instabilities.py index fc10e4d..9f85eca 100644 --- a/plot_test_secondary_instabilities.py +++ b/plot_test_secondary_instabilities.py @@ -52,7 +52,7 @@ integrator = 'classic' # first test clv -min_epsilon=10**0 +min_epsilon=10**-5 max_epsilon=101 @@ -96,8 +96,8 @@ normerror = np.memmap(savename+'/normerror_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') normerror_1st = np.memmap(savename+'/normerror_1st_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') normerrorrel = np.memmap(savename+'/normerrorrel_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') - normerrorrel_1st = np.memmap(savename+'/normerrorrel_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') - normerrorrel_2nd = np.memmap(savename+'/normerrorrel_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') + normerrorrel_1st = np.memmap(savename+'/normerrorrel_1st_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') + normerrorrel_2nd = np.memmap(savename+'/normerrorrel_2nd_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') normnonlin = np.memmap(savename+'/normnonlin_clv'+str(clv)+'.dat',mode='r',shape=(len(timeintervall),len(epsilons),len(intsteps)-1),dtype=precisioncorr,order = 'F') @@ -210,37 +210,31 @@ -# fig=plt.figure() -# for eps,epsval in enumerate(epsilons[imin:imax-1]): -# plt.plot(np.mean(np.log10(np.abs(normerrorrel[:,imin+eps,:])), axis = 0)/np.log10(epsval),np.mean(np.abs(correlationv2[:,imin+eps,:]),axis =0),label=r'$\log_{10}(\epsilon)='+str(np.log10(epsval))+r'$') -# plt.legend(loc = 'upper left',fontsize=6) -# plt.ylabel('Correlation with both orders vs rel error of both orders') -# plt.xlabel(r'$\log_{\epsilon}\left(||error||\right)$') -# plt.ylim(0,1.0) -# plt.xticks(np.arange(-8,4)) -# plt.xlim(-8,4) -# ax= plt.gca() -# ax.spines['right'].set_visible(False) -# ax.spines['top'].set_visible(False) -# fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorgrowthvscorrelation.pdf") -# fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorgrowthvscorrelation.png", dpi =400) -# plt.close(fig) + fig=plt.figure() + for ist,istval in enumerate(intsteps[1::20]): + #print(ist,istval) + plt.plot(epsilons[imin:imax],np.mean(normerrorrel_1st[:,imin:imax,ist], axis = 0),label=r'$\tau = '+str(istval*dtau)+'$') + plt.legend(loc = 'upper left',fontsize=6) + plt.title('Relative error scaling with 1st order; CLV '+str(clv),fontsize=8) + plt.xscale('log') + plt.yscale('log') + plt.grid(True) + fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_rel_1st.pdf") + fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_rel_1st.png", dpi =400) + plt.close(fig) -# fig=plt.figure() -# for eps,epsval in enumerate(epsilons[imin:imax-1]): -# plt.plot(np.mean(np.log10(np.abs(normerrorrel_1st[:,imin+eps,:])), axis = 0)/np.log10(epsval),np.mean(np.abs(correlation[:,imin+eps,:]),axis =0),label=r'$\log_{10}(\epsilon)='+str(np.log10(epsval))+r'$') -# plt.legend(loc = 'upper left',fontsize=6) -# plt.ylabel('Correlation with 1st order vs rel error of 1st order') -# plt.xlabel(r'$\log_{\epsilon}\left(||error||\right)$') -# plt.ylim(0,1.0) -# plt.xticks(np.arange(-8,4)) -# plt.xlim(-8,4) -# ax= plt.gca() -# ax.spines['right'].set_visible(False) -# ax.spines['top'].set_visible(False) -# fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorgrowthvscorrelation_1st.pdf") -# fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorgrowthvscorrelation_1st.png", dpi =400) -# plt.close(fig) + fig=plt.figure() + for ist,istval in enumerate(intsteps[1::20]): + #print(ist,istval) + plt.plot(epsilons[imin:imax],np.mean(normerrorrel[:,imin:imax,ist], axis = 0),label=r'$\tau = '+str(istval*dtau)+'$') + plt.legend(loc = 'upper left',fontsize=6) + plt.title('Relative error scaling with 1st and 2nd order; CLV '+str(clv),fontsize=8) + plt.xscale('log') + plt.yscale('log') + plt.grid(True) + fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_rel_1stand2nd.pdf") + fig.savefig(resultsfolder+"/erroranalysis/CLV_"+str(clv)+"_errorscaling_rel_1stand2nd.png", dpi =400) + plt.close(fig) # # fig=plt.figure() # for eps,epsval in enumerate(epsilons[imin:imax-1]):