在 Python 中使用“预分配单元数组”时遇到很多麻烦

Question

我最近从大量使用 MATLAB 编程转向使用 Python 编程。因此，我在运行我编写的 Python 代码时遇到了一些问题。我正在使用 numPy 和 SciPy 使用 vode 方法整合任意一组常微分方程。我已经为任意数量的“N”推广了 ODE 系统，但是遇到了预分配和使用数组的问题。

这对我来说尤其令人沮丧，因为我有两个版本的功能齐全的 MATLAB 代码，但需要将其转换为 Python 以获得优化的结果。我遇到了麻烦。尤其是以下几行：

S = np.array(np.zeros((N/2+1,1)), dtype = 'object')
KS = np.array(np.zeros((N/2+1,1)), dtype = 'object')
PS = np.array(np.zeros((N/2+1,1)), dtype = 'object')
Splot = np.array(np.zeros((N/2+1,1)), dtype = 'object')
KSplot = np.array(np.zeros((N/2+1,1)), dtype = 'object')
PSplot = np.array(np.zeros((N/2+1,1)), dtype = 'object')

代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

N = 10
K00 = np.logspace(0,3,101,10)
len1 = len(K00)
epsilon = 0.01
y0 = [0]*(3*N/2+3)
u1 = 0
u2 = 0
u3 = 0
Kplot = np.zeros((len1,1))
Pplot = np.zeros((len1,1))
S = np.array(np.zeros((N/2+1,1)), dtype = 'object')
KS = np.array(np.zeros((N/2+1,1)), dtype = 'object')
PS = np.array(np.zeros((N/2+1,1)), dtype = 'object')
Splot = np.array(np.zeros((N/2+1,1)), dtype = 'object')
KSplot = np.array(np.zeros((N/2+1,1)), dtype = 'object')
PSplot = np.array(np.zeros((N/2+1,1)), dtype = 'object')

for alpha in range(0,(N/2+1)):
    Splot[alpha] = np.zeros((len1,1))
for beta in range((N/2)+1,N+1):
    KSplot[beta-N/2-1] = np.zeros((len1,1))
for gamma in range(N+1,3*N/2+1):
    PSplot[gamma-N] = np.zeros((len1,1))

for series in range(0,len1):
    K0 = K00[series]
    Q = 10
    r1 = 0.0001
    r2 = 0.001
    a = 0.001
    d = 0.001
    k = 0.999
    S10 = 1e5
    P0 = 1

    def f(y, t):
        for alpha in range(0,(N/2+1)):
            S[alpha] = y[alpha]
        for beta in range((N/2)+1,N+1):
            KS[beta-N/2-1] = y[beta]
        for gamma in range(N+1,3*N/2+1):
            PS[gamma-N] = y[gamma]
        K = y[3*N/2+1]
        P = y[3*N/2+2]

        ydot = np.zeros((3*N/2+3,1))
        B = range((N/2)+1,N+1)
        G = range(N+1,3*N/2+1)
        runsumPS = 0
        runsum1 = 0
        runsumKS = 0 
        runsum2 = 0

        for m in range(0,N/2):
            runsumPS = runsumPS + PS[m+1]
            runsum1 = runsum1 + S[m+1]
            runsumKS = runsumKS + KS[m]
            runsum2 = runsum2 + S[m]    
            ydot[B[m]] = a*K*S[m]-(d+k+r1)*KS[m]

        for i in range(0,N/2-1):
            ydot[G[i]] = a*P*S[i+1]-(d+k+r1)*PS[i+1]

        for p in range(1,N/2):
            ydot[p] = -S[p]*(r1+a*K+a*P)+k*KS[p-1]+ \
                      d*(PS[p]+KS[p])

        ydot[0] = Q-(r1+a*K)*S[0]+d*KS[0]+k*runsumPS
        ydot[N/2] = k*KS[N/2-1]-(r2+a*P)*S[N/2]+ \
                    d*PS[N/2]
        ydot[G[N/2-1]] = a*P*S[N/2]-(d+k+r2)*PS[N/2]
        ydot[3*N/2+1] = (d+k+r1)*runsumKS-a*K*runsum2
        ydot[3*N/2+2] = (d+k+r1)*(runsumPS-PS[N/2])- \
                        a*P*runsum1+(d+k+r2)*PS[N/2]

        for j in range(0,3*N/2+3):
            return ydot[j] 

    if __name__ == '__main__':

        r = integrate.ode(f).set_integrator('vode', method='bdf')  
        t_start = 0.0
        t_final = 1e10
        delta_t = t_final/(len1-1)
        num_steps = np.floor((t_final - t_start)/delta_t) + 1

        y0[0] = S10
        for i in range(1,3*N/2+1):
            y0[i] = 0
        y0[3*N/2+1] = K0
        y0[3*N/2+2] = P0
        r.set_initial_value(y0, t_start)

        t = np.zeros((num_steps, 1))
        soln = np.array(np.zeros((num_steps, 1))*(3*N/2+3))
        t[0] = t_start
        for i in range(0,3*N/2+3):
            soln[i] = y0[i]

        k = 1
        while r.successful() and k < num_steps:
            r.integrate(r.t + delta_t)

            t[k] = r.t
            for jj in range(0,3*N/2+3):
                soln[k] = r.y[jj]
            k += 1

错误信息如下：

ValueError                                Traceback (most recent call last)
C:\Users\dis_YO_boi\Documents\Programming\Python\ArrayMod.py in <module>()
     21 
     22 for alpha in range(0,(N/2+1)):
---> 23     Splot[alpha] = np.zeros((len1,1))
     24 for beta in range((N/2)+1,N+1):
     25     KSplot[beta-N/2-1] = np.zeros((len1,1))

ValueError: could not broadcast input array from shape (101,1) into shape (1)

提前感谢您的帮助。

Answer 1

在尝试完全理解您的代码之前，让我观察一下这样的结构：

S = np.array(np.zeros((N/2+1,1)), dtype = 'object')

不好numpy。

您可能正在模仿 MATLAB 元胞数组。 Python 早在 MATLAB 之前就有“细胞”。 Python 列表可以包含各种值、字符串、数字、其他列表、数组等。

numpy 与 dtype=object 的数组只是美化的列表。如果您想要二维集合，它们可能会很方便，但与 MATLAB 一样，您无法跨此类数组的元素进行数学运算。充其量你可以迭代它们，就像你使用列表一样。

您的错误可能与此无关，但我必须深入挖掘才能确定。

---

Splot=np.array(np.zeros((4,1)),dtype=object)

产生一个(4,1) 数组，dtype=object。 np.array 尝试从输入中创建尽可能高的维度数组。

从 对于范围内的 alpha(0,(N/2+1))： Splot[alpha] = np.zeros((len1,1))

看起来您想用这些N/2+1 插槽预分配一个数组，并用一个二维数组填充每个数组。 dtype=object 有点棘手。

Splot = [np.zeros((len1,1)) for alpha in range(M)]

将生成一个包含 M 个数组的列表，每个数组的长度相同。

例如

In [67]: Splot=[np.zeros((4,1)) for alpha in range(3)]

In [68]: Splot
Out[68]: 
[array([[ 0.],
       [ 0.],
       [ 0.],
       [ 0.]]),
 array([[ 0.],
       [ 0.],
       [ 0.],
       [ 0.]]),
 array([[ 0.],
       [ 0.],
       [ 0.],
       [ 0.]])]

请注意，我将这个 2d 数组列表包装在一个数组中，我得到一个 3d 数组：

In [69]: np.array(Splot)
Out[69]: 
array([[[ 0.],
        [ 0.],
        [ 0.],
        [ 0.]],

       [[ 0.],
        [ 0.],
        [ 0.],
        [ 0.]],

       [[ 0.],
        [ 0.],
        [ 0.],
        [ 0.]]])

In [70]: _.shape
Out[70]: (3, 4, 1)

可以用这种策略制作一个数组数组。

创建一个大小合适的对象数组（np.empty 将用None 填充它）。

In [72]: Splot
Out[72]: array([0, 0, 0], dtype=object)

然后迭代以将每个0 替换为一个新对象。

In [73]: for i in range(3):
   ....:     Splot[i] = np.zeros((4,1))
   ....:     

In [74]: Splot
Out[74]: 
array([array([[ 0.],
       [ 0.],
       [ 0.],
       [ 0.]]),
       array([[ 0.],
       [ 0.],
       [ 0.],
       [ 0.]]),
       array([[ 0.],
       [ 0.],
       [ 0.],
       [ 0.]])], dtype=object)