索引錯誤 - Python，Numpy，MatLab

Question

我已經使用numpy和scipy庫將一段MatLab代碼轉換爲Python。然而，我堅持以下索引錯誤;

IndexError: index 698 is out of bounds for axis 3 with size 2 

698是時間表的大小。

它發生在代碼的最後部分，在這條線;

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

其餘的內容是爲了完整而包括在內。

以下是代碼;

import numpy as 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 into 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) 
# 
# To 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 * A 
    return y 
#======================================================================== 
def calc_expect(rhotots,Pe,nargout=1): 
    for indx in (1,2): 
     for 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 A 
#======================================================================== 


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':  
     for i in (1,(len(t)-1)): 
      if i == 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': 
     for i in (1,(len(t)-1)): 
      if i == 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) 
    for indx in (1,d1): 
     for jndx in (1,d2): 
      Aout[:,:,indx,jndx]=expm(Ain[:,:,indx,jndx]) 
    return Aout 
""" 
#======================================================================== 

#======================================================================== 


Ax=single_photon_liouvillian 
Bx=single_photon_stochastic 





toc = time.clock() 


for indx in (1,range(mt)): 
    for 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) 



for i in (2,len(t)-1): 
    #Master equation 
    rhotot=rho + dt * single_photon_liouvillian(S,L,H,rho,xi[i - 1]) 

    for indx in (1,range(mt)): 
     for 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) 




for 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) 
    for indx in (1,range(2)): 
     for 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 
    for i 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]) 
     for indx in (1,range(mt)): 
      for 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) 

任何幫助將不勝感激，因爲我一直堅持這一段時間。

謝謝！

Answer 1

的問題是，你是不是在循環定義i。 i正在被最後一次使用時遺留下來。具體來說，它使用前一個循環的最後一個值i。該值爲len(t)-1，它比長度爲2的exp_Pes的第三維的長度大得多。您需要在某處爲i定義有效值。

如果你不想遍歷i，通常你可以只在循環開始之前定義一次。但是，對於您的情況，您需要在每次循環時對其進行設置，因爲它也在此循環內定義在導致錯誤的那一行之下幾行。

但是，它可能會更清晰，只需執行類似exp_Pes[indx,jndx,j,0] = ...的操作，因爲讀取代碼的人不需要回頭查找並找到i設置的值。

一對夫婦的意見比較小的點：

np.arange不包括最後一個值。因此，在開始使用7-dt的情況下，開始時爲0，最後爲7-dt-dt，這可能不是您想要的。但無論如何，你絕對不應該在浮點數上使用np.arange（或MATLAB中的等效函數），因爲它可能存在浮點錯誤。改爲使用np.linspace。

二，對於你的_min功能，你不需要[]，你可以只做return np.amin(x), np.argmin(x)。

但是，數組通常使用方法版本更容易，方法版本是附加到數組的版本。這將是return x.min(), x.argmin()。與np.copy(rho)相同，請改爲使用rho.copy()。

但是，此功能可以進一步簡化爲lambda表達式，相當於一個Matlab匿名函數，如下所示：_min = lambda x: [x.min(), x.argmin()]（您確實需要[]）。

，您應該使用numpy版本的功能與numpy陣列。他們會更快。所以np.abs(t-t[0])比abs(t-t[0])快。與使用np.sqrt(sig1)而不是math.sqrt(sig1)一樣。

還不需要[] in returned values, so [min_dif，array_pos_start] = _min（ABS（T -t [0]））can be min_dif，array_pos_start = _min（np.abs（T-T [0]））。但是，由於您從未實際使用min_dif，因此可能只是array_pos_start = np.abs(t-t[0]).argmax()。

您不需要結束;。

在返回之前，您不需要指定給變量。例如，您可以擁有return A*B-B*A。

您似乎不必要地使用nargout參數。 Python不使用它，很大程度上是因爲你可以索引結果所以說你有一個函數foo(x)，返回最小值和最大值x。你可以只想要最大，你可以做foo(x)[1]。 Matlab不會讓你這樣做，所以它依靠nargout來解決它。