中的值错误优化曲线拟合不能通过输入numpy数组来解决

Question

我试图使用scipy.optimize.curve_fit将一些数据放入一个函数中，但是仍然会得到一个值错误，类似于this post中描述的错误，但是按照注释中的建议将所有输入输入到numpy数组中似乎并不能解决我的问题。在

以下是我的基本代码：

def RamanHamiltonian(k, omega, delta, epsilon):
    H = np.array([[(k+2.0)**2.0 - delta, omega/2.0, 0.0],
                  [omega/2.0, k**2.0-epsilon,omega/2.0],
                  [0.0,omega/2.0,(k-2.0)**2.0 + delta]])
    return H

def propagateHamiltonianTest(t, omega, delta, epsilon):  
    k = 0.0
    psi0 = np.array([0+1j*0.0, 1.0+1j*0.0, 0.0+1j*0.0])
    H = RamanHamiltonian(k, omega, delta ,epsilon)
    Energy, V = LA.eig(H)

    V = V + 1j*0.0
    Vinv = np.conjugate(np.transpose(V))

    U = np.diag(np.exp(-1j*np.array(Energy)*t))
    a  =np.dot(Vinv,psi0)
    b = np.dot(U,a)
    psi = np.dot(V,b)
    pop0 = np.absolute(psi[0])**2.0

    return pop0

popt, pcov = optimize.curve_fit(propagateHamiltonianTest,
                               np.array(tRecoils), 
                               np.array(frac0), 
                               p0=(3.0,0.05,0.03))

下面是tRecoils和frac0的值，它们的长度都是24个数组：

^{pr2}$

我得到了这个错误：

U = np.diag(np.exp(-1j*np.array(Energy)*t))

ValueError: operands could not be broadcast together with shapes (3) (24)

因此，出于某些原因，curve_nfit读取整个tRecoil数组，而不是按元素方式进行读取，我似乎无法通过更改输入格式来纠正这一点。在

Answer 1

我认为问题在于如何计算函数propagateHamiltonianTest存在歧义。它可以应用于t的单个值（这就是您对它的看法），也可以应用于整个数组，在这种情况下，numpy将整个数组传递到参数t并尝试向量化。在

事实上，当我在数组tRecoil[0]的一个元素上测试你的函数时，它可以工作，但是当你在整个数组上尝试它时，它失败了，并出现了相同的错误。所以问题与curve_fit无关。如果你可以让你的函数接受整个数组，然后返回一个数组，那么它应该与curve_fit一起工作。在

为了做到这一点，我不得不使用一个小技巧，使用scipy.linalg.block_diag。这使得我们可以创建一个更大的块对角矩阵，它基本上有每个3x3矩阵，这些矩阵是沿着对角线的每个时间步产生的。在

下面是您的propagateHamiltonianTest函数的修改版本及其注释。我已经测试过了，应该可以用了。在

from scipy.linalg import block_diag

def propagateHamiltonianTest(t, omega, delta, epsilon):  
    k = 0.0
    psi0 = np.array([0+1j*0.0, 1.0+1j*0.0, 0.0+1j*0.0])
    H = RamanHamiltonian(k, omega, delta ,epsilon)
    Energy, V = LA.eig(H)

    V = V + 1j*0.0
    Vinv = np.conjugate(np.transpose(V))

    # np.outer(t, Energy).flatten() creates a matrix for all t
    U = np.diag(np.exp(-1j*np.outer(t, Energy).flatten()))  
    a = np.dot(Vinv, psi0)
    # This repeats a so that the shape is consitent with U
    aa = block_diag(*([a]*t.size))                          
    # Have to add the transpose to make shapes match 
    b = np.dot(U, aa.T)                                     
    # Same block diagonal trick for eigenvector matrix
    VV = block_diag(*([V]*t.size))                          
    psi = np.dot(VV, b)
    pop0 = np.absolute(psi)**2.0                       
    # Since you want the first value, need to take every 3rd row 
    # and extract the values you want from the diagonal
    return np.diag(pop0[::3])

输出：

^{pr2}$