Index error - Python, Numpy, MatLab -

i have converted section of matlab code python using numpy , scipy libraries. stuck on following index error;

indexerror: index 698 out of bounds axis 3 size 2 

698 size of time list.

it occurs in last section of code, on line;

exp_pes[indx,jndx,j,i]=np.trace(rhotemp[:,:,indx,jndx] * pe) 

the rest included completeness.

the following code;

import numpy np import math import time #from scipy.linalg import expm tic = time.clock()  dt = 1e-2;                      # time step t = np.arange(0, 7-dt, dt)         # t   gam=1 mt=2 nt=2 delta=1     ensemble_size=6  rkflag=0   itype = 'em'     #======================================================================== def _min(x):     return [np.amin(x), np.argmin(x)] #======================================================================== ############### #wavepacket envelope #rising exponential  sig1 = 1; xi_t = np.zeros(len(t)); [min_dif, array_pos_start] = _min(abs(t -t[0] )); [min_dif, array_pos_stop ] = _min(abs(t -t[-1]/2.0)); step = t[1]-t[0] p_len = np.arange(0, t[array_pos_stop] - t[array_pos_start] + step, step) xi_t[array_pos_start:array_pos_stop+1] = math.sqrt(sig1) * np.exp((sig1/2.0)*p_len); norm = np.trapz(t, xi_t * np.conj(xi_t)); xi = xi_t/np.sqrt(abs(norm));   #pauli matrices #======================== pe=np.array([[1,0],[0,0]]) sm=np.array([[0,0],[1,0]]) sz=np.array([[1,0],[0,- 1]]) sy=np.array([[0,- 1j],[1j,0]]) sx=np.array([[0,1],[1,0]])   #s,l,h coefficients #========================= s=np.eye(2) l=np.sqrt(gam) * sm h=np.eye(2) #=========================     psi=np.array([[0],[1]]) rhoi=psi * psi.t  #rho=np.zeros(2,2,2,2) rho=np.zeros((2,2,2,2)) rhotot=np.copy(rho) rhoblank=np.copy(rho) rhos=np.copy(rho) rhotots=np.copy(rho) rhon=np.copy(rho) rhototn=np.copy(rho) rhov=np.copy(rhoi)  exp_pes=np.zeros((2,2,2,2)) exp_pes2=np.zeros((2,2,2,2)) #initial conditions rho #==================  rho[:,:,0,0]=rhoi #rho[:,:,2,2]=rhoi rhosi=np.copy(rho)    num_bad=0 avg_val=0 num_badrk=0 avg_valrk=0  #np.zeros(len(t));  exp_sz=np.zeros((2,2,len(t))) exp_sy=np.copy(exp_sz) exp_sx=np.copy(exp_sz) exp_pe=np.copy(exp_sz)    ##########################3 #functions      #======================================================================== #d[x]rho = x*rho* ctranspose(x) -0.5*(ctranspose(x)*x*rho + rho*ctranspose(x)*x) # # call write curlyd(x,rho)  def curlyd(x,rho,nargout=1):     y=x * rho * x.conj().t - 0.5 * (x.conj().t * x * rho + rho * x.conj().t * x)     return y  #======================================================================== def curlyc(a,b,nargout=1):     y=a * b - b *     return y #======================================================================== def calc_expect(rhotots,pe,nargout=1):     indx in (1,2):         jndx in (1,2):             exp_pes[indx,jndx]=np.trace(rhotots[:,:,indx,jndx] * pe)             exp_pes2[indx,jndx]=np.trace(rhotots[:,:,indx,jndx] * pe) ** 2     return exp_pes,exp_pes2 #========================================================================  #========================================================================     #========================================================================  #========================================================================          #========================================================================  def single_photon_liouvillian(s,l,h,rho,xi,nargout=1):     rhotot[:,:,1,1]=curlyd(l,rho[:,:,1,1]) + xi * curlyc(s * rho[:,:,0,1],l.t) + xi.t * curlyc(l,rho[:,:,1,0] * s.t) + xi.t * xi * (s * rho[:,:,0,0] * s.t - rho[:,:,0,0])     rhotot[:,:,1,0]=curlyd(l,rho[:,:,1,0]) + xi * curlyc(s * rho[:,:,0,0],l.t)     rhotot[:,:,0,1]=curlyd(l,rho[:,:,0,1]) + xi.t * curlyc(l,rho[:,:,0,0] * s.t)     rhotot[:,:,0,0]=curlyd(l,rho[:,:,0,0])     a=np.copy(rhotot)     return #========================================================================   def single_photon_stochastic(s,l,h,rho,xi,nargout=1):     k=np.trace((l + l.t) * rho[:,:,1,1]) + np.trace(s * rho[:,:,0,1]) * xi + np.trace(s.t * rho[:,:,1,0]) * xi.t     rhotot[:,:,1,1]=(l * rho[:,:,1,1] + rho[:,:,1,1] * l.t + s * rho[:,:,0,1] * xi + rho[:,:,1,0] * s.t * xi.t - k * rho[:,:,1,1])     rhotot[:,:,1,0]=(l * rho[:,:,1,0] + rho[:,:,1,0] * l.t + s * rho[:,:,0,0] * xi - k * rho[:,:,1,0])     rhotot[:,:,0,1]=(l * rho[:,:,0,1] + rho[:,:,0,1] * l.t + rho[:,:,0,0] * s.t * xi.t - k * rho[:,:,0,1])     rhotot[:,:,0,0]=(l * rho[:,:,0,0] + rho[:,:,0,0] * l.t - k * rho[:,:,0,0])     b=np.copy(rhotot)     return b  #======================================================================== def sde_int_io2_rk(a,b,s,l,h,yn,xi,dt,dw,nargout=1):     gn=yn + a[s,l,h,yn,xi] * dt + b[s,l,h,yn,xi] * dw     gnp=yn + a[s,l,h,yn,xi] * dt + b[s,l,h,yn,xi] * np.sqrt(dt)     gnm=yn + a[s,l,h,yn,xi] * dt - b[s,l,h,yn,xi] * np.sqrt(dt)     ynp1=yn + 0.5 * (a[s,l,h,gn,xi] + a[s,l,h,yn,xi]) * dt + 0.25 * (b[s,l,h,gnp,xi] + b[s,l,h,gnm,xi] + 2 * b[s,l,h,yn,xi]) * dw + 0.25 * (b[s,l,h,gnp,xi] - b[s,l,h,gnm,xi]) * (dw ** 2 - dt) * (dt) ** (- 0.5)     return ynp1 #========================================================================    #======================================================================== def sde_int_photon(itype,rhos,s,l,h,pe,xi,t,nargout=1):     dt=t[1] - t[0]     rhoblank=np.zeros(len(rhos))     ax=single_photon_liouvillian     bx=single_photon_stochastic     #if strcmp(itype,'em'):     if itype == 'em':             in (1,(len(t)-1)):             if == 1:                 exp_pes[:,:,i],exp_pes2[:,:,i]=calc_expect(rhos,pe,nargout=2)                 continue             dw=np.sqrt(dt) * np.randn             rhotots=rhos + dt * single_photon_liouvillian(s,l,h,rhos,xi[i]) + dw * single_photon_stochastic(s,l,h,rhos,xi[i])             exp_pes[:,:,i],exp_pes2[:,:,i]=calc_expect(rhotots,pe,nargout=2)             rhos=np.copy(rhotots)             rhotots=np.copy(rhoblank)     if itype == 'rk':          in (1,(len(t)-1)):             if == 1:                 exp_pes[:,:,i],exp_pes2[:,:,i]=calc_expect(rhos,pe,nargout=2)                 continue             dw=np.sqrt(dt) * np.randn             rhotots=sde_int_io2_rk(ax,bx,s,l,h,rhos,xi[i],dt,dw)             exp_pes[:,:,i],exp_pes2[:,:,i]=calc_expect(rhotots,pe,nargout=2)             rhos=np.copy(rhotots)             rhotots=np.copy(rhoblank)     return exp_pes,exp_pes2 #========================================================================  """ def md_expm(ain,nargout=1):     aout=np.zeros(len(ain))     r,c,d1,d2=len(ain,nargout=4)     indx in (1,d1):         jndx in (1,d2):             aout[:,:,indx,jndx]=expm(ain[:,:,indx,jndx])     return aout """ #========================================================================  #========================================================================   ax=single_photon_liouvillian bx=single_photon_stochastic      toc = time.clock()   indx in (1,range(mt)):     jndx in (1,range(nt)):         exp_sz[indx,jndx,1]=np.trace(rho[:,:,indx,jndx] * sz)         exp_sy[indx,jndx,1]=np.trace(rho[:,:,indx,jndx] * sy)         exp_sx[indx,jndx,1]=np.trace(rho[:,:,indx,jndx] * sx)         exp_pe[indx,jndx,1]=np.trace(rho[:,:,indx,jndx] * pe)    in (2,len(t)-1):     #master equation     rhotot=rho + dt * single_photon_liouvillian(s,l,h,rho,xi[i - 1])      indx in (1,range(mt)):         jndx in (1,range(nt)):             exp_sz[indx,jndx,i]=np.trace(rhotot[:,:,indx,jndx] * sz)             exp_sy[indx,jndx,i]=np.trace(rhotot[:,:,indx,jndx] * sy)             exp_sx[indx,jndx,i]=np.trace(rhotot[:,:,indx,jndx] * sx)             exp_pe[indx,jndx,i]=np.trace(rhotot[:,:,indx,jndx] * pe)     rho=np.copy(rhotot)     rhotot=np.copy(rhoblank)     j in (1,range(ensemble_size)):     psi1=np.array([[0],[1]])     rho1=psi1 * psi1.t     rhotemp = np.zeros((2,2,2,2))     rhotemp[:,:,0,0]=rho1     rhotemp[:,:,1,1]=rho1     rhos=np.copy(rhotemp)     indx in (1,range(2)):         jndx in (1,range(2)):             exp_pes[indx,jndx,j,i]=np.trace(rhotemp[:,:,indx,jndx] * pe)             exp_pes2[indx,jndx,j,i]=np.trace(rhotemp[:,:,indx,jndx] * pe) ** 2     in (2,(len(t)-1)):         dw=np.sqrt(dt) * np.random.randn()         rhotots=rhos + dt * single_photon_liouvillian(s,l,h,rhos,xi[i - 1]) + dw * single_photon_stochastic(s,l,h,rhos,xi[i - 1])         indx in (1,range(mt)):             jndx in (1,range(nt)):                 exp_pes[indx,jndx,j,i]=np.trace(rhotots[:,:,indx,jndx] * pe)                 exp_pes2[indx,jndx,j,i]=np.trace(rhotots[:,:,indx,jndx] * pe) ** 2         rhos=np.copy(rhotots)         rhotots=np.copy(rhoblank)     irow=np.where(np.squeeze(exp_pes[2,2,j,:]) > 1)     val=np.squeeze(exp_pes[2,2,j,irow])     if irow:         num_bad=num_bad + 1         avg_val=avg_val + max(val) 

any appreciated have been stuck on while.

thanks!

the problem not defining i in loop. i being left on last time used. specifically, using last value of i previous loop. value len(t)-1, larger length of 3rd dimension of exp_pes has length of 2. need define valid value i somewhere.

if don't want loop on i, define once before loop starts. in case, however, need set every time loop since being defined inside loop few lines below 1 causing error.

however, clearer exp_pes[indx,jndx,j,0] = ..., since people reading code won't need , find value i set to.

a couple more minor points of advice:

np.arange not include last value. case use 7-dt @ beginning, ending start being 0 , end being 7-dt-dt, not want. anyway, should never use np.arange (or equivalent in matlab) floats, since can have floating-point errors. use np.linspace instead.

second, _min function, don't need [], can return np.amin(x), np.argmin(x).

however, easier arrays use method version, version attached array. return x.min(), x.argmin(). same np.copy(rho), use rho.copy() instead.

however, function further simplified lambda expression, equivalent of matlab anonymous function, so: _min = lambda x: [x.min(), x.argmin()] (you need [] there).

you should use numpy version of functions numpy arrays. faster. np.abs(t-t[0]) faster abs(t-t[0]). same using np.sqrt(sig1) instead of math.sqrt(sig1).

you don't need [] in returned values, so[min_dif, array_pos_start] = _min(abs(t -t[0] ))can bemin_dif, array_pos_start = _min(np.abs(t-t[0])). however, since never using min_dif, array_pos_start = np.abs(t-t[0]).argmax().

you don't need end lines ;.

you don't need assign variable before returning. can have, example, return a*b-b*a.

you seem unnecessarily use nargout argument. python doesn't use this, in large part because can index results have function foo(x) returns min , max of x. can want max, can foo(x)[1]. matlab doesn't let that, relies on nargout work around it.