3

我有一个m x n数组：a，其中整数m > 1E6和n <= 5。

我有函数F和G，它们的组成如下：F ( u , G ( u , t))。u是一个1 x n数组，t 是一个标量，F和G返回1 x n数组。

我需要评估F中的每row一个，并将先前评估的行用作下一次评估的u -array。我需要做出这样的评价。am

这必须非常快。我之前对scitools.std StringFunction整个数组的评估印象深刻，但是这个问题需要使用之前计算的数组作为计算下一个数组的参数。我不知道 StringFunction 是否可以做到这一点。

例如：

a = zeros((1000000, 4))
a[0] = asarray([1.,69.,3.,4.1])

# A is a float defined elsewhere, h is a function which accepts a float as its argument and returns an arbitrary float. h is defined elsewhere.

def G(u, t):
  return asarray([u[0], u[1]*A, cos(u[2]), t*h(u[3])])

def F(u, t):
  return u + G(u, t)


dt = 1E-6

for i in range(1, 1000000):
  a[i] = F(a[i-1], i*dt)
  i += 1

上面代码的问题是它慢得要命。我需要在 numpy 毫秒内完成这些计算。

我怎样才能做我想做的事？

谢谢你的时间。

亲切的问候，

马吕斯

4

2 回答 2

1

这种事情在 numpy 中很难做到。如果我们按列查看，我们会看到一些更简单的解决方案。

a[:,0]很容易：

col0 = np.ones((1000))*2
col0[0] = 1                  #Or whatever start value.
np.cumprod(col0, out=col0)

np.allclose(col0, a[:1000,0])
True

如前所述，这将很快溢出。a[:,1]可以沿着同样的路线做很多事情。

我不相信有办法单独在 numpy 中快速完成接下来的两列。我们可以为此转向 numba：

from numba import auotojit

def python_loop(start, count):
     out = np.zeros((count), dtype=np.double)
     out[0] = start
     for x in xrange(count-1):
         out[x+1] = out[x] + np.cos(out[x+1])
     return out

numba_loop = autojit(python_loop)

np.allclose(numba_loop(3,1000),a[:1000,2])
True

%timeit python_loop(3,1000000)
1 loops, best of 3: 4.14 s per loop

%timeit numba_loop(3,1000000)
1 loops, best of 3: 42.5 ms per loop

尽管值得指出的是，它收敛pi/2得非常快，并且对于任何起始值，计算超过约 20 个值的递归几乎没有意义。这将返回与双点精度完全相同的答案——我没有费心找到截止值，但它远小于 50：

%timeit tmp = np.empty((1000000)); 
        tmp[:50] = numba_loop(3,50);
        tmp[50:] = np.pi/2
100 loops, best of 3: 2.25 ms per loop

您可以对第四列执行类似的操作。当然，您可以autojit使用所有功能，但这会根据 numba 的使用情况为您提供几种不同的选择：

  1. 对前两列使用 cumprod
  2. 对第 3 列（和可能的第 4 列）使用近似值，其中仅计算前几次迭代
  3. 使用 numba 实现第 3 列和第 4 列autojit
  4. 将所有内容包装在 autojit 循环中（最佳选择）
  5. 您呈现超过 ~200 的所有行的方式将是np.inf或np.pi/2. 利用这一点。
于 2013-12-17T22:17:28.147 回答
0

稍微快一点。您的第一列基本上是 2^n。为 n 计算 2^n 高达 1000000 会溢出。第二列更糟。

def calc(arr, t0=1E-6):
    u = arr[0]
    dt = 1E-6
    h = lambda x: np.random.random(1)*50.0

    def firstColGen(uStart):
        u = uStart
        while True:
            u += u
            yield u

    def secondColGen(uStart, A):
        u = uStart
        while True:
            u += u*A
            yield u

    def thirdColGen(uStart):
        u = uStart
        while True:
            u += np.cos(u)
            yield u

    def fourthColGen(uStart, h, t0, dt):
        u = uStart
        t = t0
        while True:
            u += h(u) * dt
            t += dt
            yield u

    first = firstColGen(u[0])
    second = secondColGen(u[1], A)
    third = thirdColGen(u[2])
    fourth = fourthColGen(u[3], h, t0, dt)

    for i in xrange(1, len(arr)):
        arr[i] = [first.next(), second.next(), third.next(), fourth.next()]
于 2013-12-17T19:31:35.283 回答