optimize.curve_fit 中的值错误无法通过输入 numpy 数组解决
Value error in optimize.curve_fit not solved by inputing numpy arrays
我正在尝试使用 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 的数组:
tRecoils = array([ 2.88597836, 1.15439135, 1.73158702, 2.19334356, 0.23087827,
2.07790442, 0.11543913, 1.50070875, 2.77053923, 1.61614788,
1.03895221, 2.42422183, 0.92351308, 0.80807394, 0.3463174 ,
0.57719567, 2.6551001 , 0.46175654, 1.84702615, 1.38526961,
2.53966096, 0.69263481, 2.30878269, 1.96246529])
frac0 = array([ 0.15761062, 0.02044625, 0.17275937, 0.02236243, 0.07388558,
0.00967176, 0.01886309, 0.20412516, 0.21667489, 0.21783697,
0.00173812, 0.14038657, 0.03145599, 0.08644404, 0.13153078,
0.18794377, 0.2139092 , 0.17141201, 0.13021916, 0.12671806,
0.21090369, 0.1611094 , 0.08732627, 0.05764911])
我收到此错误:
U = np.diag(np.exp(-1j*np.array(Energy)*t))
ValueError: operands could not be broadcast together with shapes (3)
(24)
所以出于某种原因 curve_fit 读取整个 tRecoil
数组而不是按元素读取,我似乎无法通过更改输入格式来解决这个问题。
我认为问题在于您的函数 propagateHamiltonianTest
的计算方式存在歧义。它可以应用于 t
的单个值(这就是您的想法),也可以应用于整个数组,在这种情况下,numpy 将整个数组传递给参数 t
并尝试向量化。
事实上,例如,当我在数组 tRecoil[0]
的单个元素上测试您的函数时,它可以工作,但是当您在整个数组上尝试它时,它会失败并出现相同的错误。所以问题与curve_fit
无关。如果你可以让你的函数将整个数组输入,return 一个数组输出,那么它应该与 curve_fit
.
一起工作
为了做到这一点,我不得不使用一些 hacky 技巧,使用 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])
输出:
print('popt: ', popt, '\n\npcov: ', pcov)
popt: [ 2.59126808 -0.13125704 0.52597681]
pcov: [[ 0.00114962 0.00072806 -0.0005113 ]
[ 0.00072806 0.00076575 -0.00020831]
[-0.0005113 -0.00020831 0.00063511]]
我正在尝试使用 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 的数组:
tRecoils = array([ 2.88597836, 1.15439135, 1.73158702, 2.19334356, 0.23087827,
2.07790442, 0.11543913, 1.50070875, 2.77053923, 1.61614788,
1.03895221, 2.42422183, 0.92351308, 0.80807394, 0.3463174 ,
0.57719567, 2.6551001 , 0.46175654, 1.84702615, 1.38526961,
2.53966096, 0.69263481, 2.30878269, 1.96246529])
frac0 = array([ 0.15761062, 0.02044625, 0.17275937, 0.02236243, 0.07388558,
0.00967176, 0.01886309, 0.20412516, 0.21667489, 0.21783697,
0.00173812, 0.14038657, 0.03145599, 0.08644404, 0.13153078,
0.18794377, 0.2139092 , 0.17141201, 0.13021916, 0.12671806,
0.21090369, 0.1611094 , 0.08732627, 0.05764911])
我收到此错误:
U = np.diag(np.exp(-1j*np.array(Energy)*t))
ValueError: operands could not be broadcast together with shapes (3) (24)
所以出于某种原因 curve_fit 读取整个 tRecoil
数组而不是按元素读取,我似乎无法通过更改输入格式来解决这个问题。
我认为问题在于您的函数 propagateHamiltonianTest
的计算方式存在歧义。它可以应用于 t
的单个值(这就是您的想法),也可以应用于整个数组,在这种情况下,numpy 将整个数组传递给参数 t
并尝试向量化。
事实上,例如,当我在数组 tRecoil[0]
的单个元素上测试您的函数时,它可以工作,但是当您在整个数组上尝试它时,它会失败并出现相同的错误。所以问题与curve_fit
无关。如果你可以让你的函数将整个数组输入,return 一个数组输出,那么它应该与 curve_fit
.
为了做到这一点,我不得不使用一些 hacky 技巧,使用 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])
输出:
print('popt: ', popt, '\n\npcov: ', pcov)
popt: [ 2.59126808 -0.13125704 0.52597681]
pcov: [[ 0.00114962 0.00072806 -0.0005113 ]
[ 0.00072806 0.00076575 -0.00020831]
[-0.0005113 -0.00020831 0.00063511]]