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

Question

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

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

мне нужно оценить каждую из rowa в F, и использовать предварительно оцененную строку в качестве у -array для следующей оценки. Мне нужно сделать 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 

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

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

Благодарим вас за наше время.

С наилучшими пожеланиями,

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] можно сделать много по тем же линиям.

Я не верю, что есть способ сделать следующие две колонки внутри 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. Реализовать столбцы 3 и 4 в Numba с помощью autojit
4. Wrap все внутри петли autojit (лучший вариант)
5. так, как вы представили это все строки мимо ~ 200 будет либо b e 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()]