具有在Python'preallocating电池阵列“的麻烦 [英] Having much trouble with 'preallocating cell arrays' in Python

问题描述

我最近从在Python在MATLAB大量编程编程切换。因此，我有跑，我已经写了Python code的一些问题。我使用numpy的和SciPy的整合常微分方程任意一组，用方法的 VODE 。我曾概括为常微分方程他们中的任何数量的'N'的系统，但碰到有preallocating和使用数组的问题。

I have recently switched from programming heavily in MATLAB to programming in Python. Hence I am having some issues running the Python code that I have written. I am using numPy and SciPy to integrate an arbitrary set of ordinary differential equations, with method vode. I have generalized the system of ODEs for any number 'N' of them, but run into problems with preallocating and using arrays.

因为我有两个版本，全功能的MATLAB code，但需要将其转换到Python优化的结果，是对我特别沮丧。而且我有麻烦。尤其是以下行：

It is especially frustrating for me because I have TWO versions of fully functional MATLAB code, but need to convert it to Python for optimized results. And I am having trouble. Especially with the following lines:

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')

下面是code：

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

该错误消息如下：

The error message is below:

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)

感谢您事先的任何帮助。

Thanks in advance for any help.

推荐答案

想完全理解你的code之前，让我观察到结构，如：

Before trying to fully understand your code, let me observe that constructs like:

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

不是很好 numpy的。

您可能正在模仿MATLAB电池阵列。蟒蛇有MATLAB过不了多久'细胞'。 Python列表可以容纳多种多样的设定值，字符串，数字，其他列表，数组的等。

You are probably imitating the MATLAB cell arrays. Python had 'cells' long before MATLAB. Python lists can hold are diverse set of values, strings, numbers, other lists, arrays, etc.

numpy的数组只是荣耀列表。他们可能是方便，如果你想2D的集合，但与MATLAB你不能碰到这样数组的元素做数学。在最好的情况你可以对他们进行迭代，就像你用列表。

numpy arrays with dtype=object are just glorified lists. They may be convenient if you want 2d collections, but as with MATLAB you can't do math across elements of such arrays. At best you can iterate over them, just as you would with lists.

这可能是因为你的错误有没有关系，但我得掏一点是肯定的。

It's possible that your error has nothing to do with that, but I'll have to dig a bit to be sure.

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

生成一个（4,1）数组，DTYPE =对象。 np.array 尝试创建尽可能高的维数组，因为它可以从输入端。

produces a (4,1) array, with dtype=object. np.array tries to create as high a dimensional array as it can from the inputs.

从
    有效范围内的α（0，（N / 2 + 1））：
        Splotα= np.zeros（（len1,1））

From for alpha in range(0,(N/2+1)): Splot[alpha] = np.zeros((len1,1))

它看起来像你想preallocate的这些 N / 2 + 1 槽的阵列，并用二维数组填充每个。这是一个有点棘手与 DTYPE =对象。

it looks like you wanted to preallocate an an array with these N/2+1 slots, and fill each with a 2d array. It's a little tricky with dtype=object.

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

会产生具有M阵列的列表，每一个具有相同的长度。

would produce a list with M arrays, each of the same length.

例如

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.]])]

请注意，我换数组中的二维数组的这个名单，我得到一个三维数组：

Note that I wrap this list of 2d arrays in an array, I get a 3d array:

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)

有可能使数组的数组这个策略

It is possible to make an array of arrays with this strategy.

请大小合适的对象数组（ np.empty 将无填充）。

Make an object array of the right size (np.empty will fill it with 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)

这篇关于具有在Python'preallocating电池阵列“的麻烦的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！