numpy: вычисление функции в матрице, использование предыдущего массива в качестве аргумента при вычислении следующего

Question

У меня есть m x n массив: aгде целые числа m > 1E6, а также n <= 5,

У меня есть функции F и G, которые составлены так: F(u, G (u, t)). ты это 1 x n массив, т является скалярным, и F и G возвращает 1 x n массивы.

Мне нужно оценить каждый row из a в F, и использовать ранее оцененную строку в качестве u- массива для следующей оценки. Мне нужно сделать m такие оценки.

Это должно быть очень быстро. Я был ранее впечатлен scitools.stdStringFunction оценка для всего массива, но эта проблема требует использования ранее вычисленного массива в качестве аргумента при расчете следующего. Я не знаю, может ли 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

Проблема с приведенным выше кодом заключается в том, что он медленный, как ад. Мне нужно, чтобы эти расчеты были сделаны за несколько миллисекунд.

Как я могу делать то, что я хочу?

Спасибо за наше время.

С уважением,

Marius

Answer 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] можно сделать много в том же духе.

Я не верю, что есть способ быстро выполнить следующие две колонки в одиночестве. Мы можем обратиться к 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. Реализуйте столбцы 3 и 4 в numba, используя autojit
4. Оберните все внутри цикла autojit (лучший вариант)
5. То, как вы представили это, все строки после ~200 будут либо np.inf или же np.pi/2, Используйте это.

Answer 2

Чуть быстрее Ваш первый столбец в основном 2^n. Вычисление 2 ^ n для 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()]