Question

T(i) = Tm(i) + (T(i-1)-Tm(i))**(-tau(i))

Tm и tau - это векторы NumPy той же длины, которые были рассчитаны ранее, и необходимо создать новый вектор T.i включено только для указания индекса элемента для того, что требуется.

Необходим ли цикл for для этого случая?

Andrew Jaffe · Answer 1 · 10 декабря 2010

Вы можете подумать, что это будет работать:

import numpy as np
n = len(Tm)
t = np.empty(n)

t[0] = 0  # or whatever the initial condition is 
t[1:] = Tm[1:] + (t[0:n-1] - Tm[1:])**(-tau[1:])

но это не так: на самом деле вы не можете сделать рекурсию в numpy таким образом (поскольку numpy вычисляет всю RHS и затем назначает ее для LHS).

Итак, если вы не можете придумать нерекурсивную версию этой формулы, вы застряли с явным циклом:

tt = np.empty(n)
tt[0] = 0.
for i in range(1,n):
    tt[i] = Tm[i] + (tt[i-1] - Tm[i])**(-tau[i])

Sven Marnach · Answer 2 · 10 декабря 2010

В некоторых случаях возможна такая рекурсия, а именно, если для формулы рекурсии есть Ufunc NumPy, например,

c = numpy.arange(10.)
numpy.add(c[:-1], c[1:], c[1:])

Это вычисляет накопленные суммы на cиспользуя выходной параметр numpy.add.Его нельзя записать как

c[1:] = c[:-1] + c[1:]

, потому что теперь результатом добавления является временный объект, который копируется в c[1:] после завершения вычислений.

Наиболее естественная вещь дляпопробуйте сейчас определить ваш собственный ufunc:

def f(T, Tm, tau):
    return Tm + (T - Tm)**(-tau)
uf = numpy.frompyfunc(f, 3, 1)

Но по не зависящим от меня причинам следующее не работает:

uf(T[:-1], Tm[1:], tau[1:], T[1:])

Видимо, результат напрямую не записывается T[1:], но хранится во временном хранилище и копируется после завершения операции.Даже если бы это сработало, я бы не ожидал от этого ускорения по сравнению с обычным циклом, поскольку для каждой записи необходимо вызывать функцию Python.

Если вы действительно хотите избежать цикла Python, вывероятно, нужно перейти на Cython или ctypes.

keiv.fly · Answer 3 · 18 мая 2016

Я выполнил некоторые тесты, и в 2018 году, используя Numba - это первый вариант, который люди должны пытаться ускорить рекурсивными функциями в Numpy (предложение Aronstef).Numba уже предустановлена в пакете Anaconda и имеет один из самых быстрых времен (примерно в 20 раз быстрее, чем любой Python).В 2018 году Python поддерживает аннотации @numba без дополнительных шагов (как минимум версии 3.6 и 3.7).Вот два теста: один выполнен 2018-10-20, а другой 2016-05-18.

И, как упоминал Джаффе, в 2018 году по-прежнему невозможно векторизовать рекурсивные функции.Я проверил векторизацию с помощью Aronstef, и она НЕ работает.

Тесты, отсортированные по времени выполнения:

-----------------------------------
|Variant        |2018-10 |2016-05 |
-----------------------------------
|Pure C         |   na   | 2.75 ms|
|C extension    |   na   | 6.22 ms|
|Cython float32 | 1.01 ms|   na   |
|Cython float64 | 1.05 ms| 6.26 ms|
|Fortran f2py   |   na   | 6.78 ms|
|Numba float32  | 2.81 ms|   na   |
|(Aronstef)     |        |        |
|Numba float64  | 5.28 ms|   na   |
|Append to list |48.2  ms|91.0  ms|
|Using a.item() |58.3  ms|74.4  ms|
|np.fromiter()  |60.0  ms|78.1  ms|
|Loop over Numpy|71.9  ms|87.9  ms|
|(Jaffe)        |        |        |
|Loop over Numpy|74.4  ms|   na   |
|(Aronstef)     |        |        |
-----------------------------------

Соответствующий код предоставляется в конце ответа.

Я не проверял Pure C в 2018 году, но я полагаю, что он все еще самый быстрый на основе предыдущего теста.

Я также не проверял расширение C в 2018 году, и я думаю, что оно почти такое же время, как на основе Cython.на предыдущем бенчмарке.

Fortran было очень сложно отлаживать и компилировать, поэтому я не проверял версию f2py в 2018 году. И все равно это было хуже, чем Cython.

У меня есть следующееНастройка в 2018 году:

Processor: Intel i7-7500U 2.7GHz
Versions:
Python:  3.7.0
Numba:  0.39.0
Cython: 0.28.5
Numpy:  1.15.1

Рекомендуемый Numba код с использованием float32 (от Aronstef):

@numba.jit("float32[:](float32[:], float32[:])", nopython=False, nogil=True)
def calc_py_jit32(Tm_, tau_):
    tt = np.empty(len(Tm_),dtype="float32")
    tt[0] = Tm_[0]
    for i in range(1, len(Tm_)):
        tt[i] = Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i])
    return tt[1:]

Все остальные коды:

Создание данных (например, комментарий Aronstef + Mike T):

np.random.seed(0)
n = 100000
Tm = np.cumsum(np.random.uniform(0.1, 1, size=n).astype('float64'))
tau = np.random.uniform(-1, 0, size=n).astype('float64')
ar = np.column_stack([Tm,tau])
Tm32 = Tm.astype('float32')
tau32 = tau.astype('float32')
Tm_l = list(Tm)
tau_l = list(tau)

Код в 2016 году немного отличался, поскольку я использовал функцию abs () для предотвращения nans, а не вариант Mike T. В 2018 году функцияточно так же, как OP (OrigИнал Плакат) писал.

Cython float32 с использованием магии Jupyter %%.Эту функцию можно использовать непосредственно в Python.Cython нужен компилятор C ++, в котором был скомпилирован Python.Установка правильной версии компилятора Visual C ++ (для Windows) может быть проблематичной:

%%cython

import cython
import numpy as np
cimport numpy as np
from numpy cimport ndarray

cdef extern from "math.h":
    np.float32_t exp(np.float32_t m)

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.infer_types(True)
@cython.initializedcheck(False)

def cy_loop32(np.float32_t[:] Tm,np.float32_t[:] tau,int alen):
    cdef np.float32_t[:] T=np.empty(alen, dtype=np.float32)
    cdef int i
    T[0]=0.0
    for i in range(1,alen):
        T[i] = Tm[i] + (T[i-1] - Tm[i])**(-tau[i])
    return T

Cython float64 с использованием Jupyter %% magic.Эту функцию можно использовать непосредственно в Python:

%%cython

cdef extern from "math.h":
    double exp(double m)
import cython
import numpy as np
cimport numpy as np
from numpy cimport ndarray

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.infer_types(True)
@cython.initializedcheck(False)

def cy_loop(double[:] Tm,double[:] tau,int alen):
    cdef double[:] T=np.empty(alen)
    cdef int i
    T[0]=0.0
    for i in range(1,alen):
        T[i] = Tm[i] + (T[i-1] - Tm[i])**(-tau[i])
    return T

Numba float64:

@numba.jit("float64[:](float64[:], float64[:])", nopython=False, nogil=True)
def calc_py_jit(Tm_, tau_):
    tt = np.empty(len(Tm_),dtype="float64")
    tt[0] = Tm_[0]
    for i in range(1, len(Tm_)):
        tt[i] = Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i])
    return tt[1:]

Добавить в список .Самое быстрое не скомпилированное решение:

def rec_py_loop(Tm,tau,alen):
     T = [Tm[0]]
     for i in range(1,alen):
        T.append(Tm[i] - (T[i-1] + Tm[i])**(-tau[i]))
     return np.array(T)

Использование a.item ():

def rec_numpy_loop_item(Tm_,tau_):
    n_ = len(Tm_)
    tt=np.empty(n_)
    Ti=tt.item
    Tis=tt.itemset
    Tmi=Tm_.item
    taui=tau_.item
    Tis(0,Tm_[0])
    for i in range(1,n_):
        Tis(i,Tmi(i) - (Ti(i-1) + Tmi(i))**(-taui(i)))
    return tt[1:]

np.fromiter ():

def it(Tm,tau):
    T=Tm[0]
    i=0
    while True:
        yield T
        i+=1
        T=Tm[i] - (T + Tm[i])**(-tau[i])

def rec_numpy_iter(Tm,tau,alen):
    return np.fromiter(it(Tm,tau), np.float64, alen)[1:]

Зацикливание на Numpy (на основе идеи Джаффа):

def rec_numpy_loop(Tm,tau,alen):
    tt=np.empty(alen)
    tt[0]=Tm[0]
    for i in range(1,alen):
        tt[i] = Tm[i] - (tt[i-1] + Tm[i])**(-tau[i])
    return tt[1:]

Зацикливание на Numpy (код Аронстефа). На моемкомпьютер float64 является типом по умолчанию для np.empty.

def calc_py(Tm_, tau_):
    tt = np.empty(len(Tm_),dtype="float64")
    tt[0] = Tm_[0]
    for i in range(1, len(Tm_)):
        tt[i] = (Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i]))
    return tt[1:]

Pure C без использования Python вообще.Версия от 2016 года (с функцией fabs ()):

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <windows.h>
#include <sys\timeb.h> 

double randn() {
    double u = rand();
    if (u > 0.5) {
        return sqrt(-1.57079632679*log(1.0 - pow(2.0 * u - 1, 2)));
    }
    else {
        return -sqrt(-1.57079632679*log(1.0 - pow(1 - 2.0 * u,2)));
    }
}
void rec_pure_c(double *Tm, double *tau, int alen, double *T)
{

    for (int i = 1; i < alen; i++)
    {
        T[i] = Tm[i] + pow(fabs(T[i - 1] - Tm[i]), (-tau[i]));
    }
}

int main() {
    int N = 100000;
    double *Tm= calloc(N, sizeof *Tm);
    double *tau = calloc(N, sizeof *tau);
    double *T = calloc(N, sizeof *T);
    double time = 0;
    double sumtime = 0;
    for (int i = 0; i < N; i++)
    {
        Tm[i] = randn();
        tau[i] = randn();
    }

    LARGE_INTEGER StartingTime, EndingTime, ElapsedMicroseconds;
    LARGE_INTEGER Frequency;
    for (int j = 0; j < 1000; j++)
    {
        for (int i = 0; i < 3; i++)
        {
            QueryPerformanceFrequency(&Frequency);
            QueryPerformanceCounter(&StartingTime);

            rec_pure_c(Tm, tau, N, T);

            QueryPerformanceCounter(&EndingTime);
            ElapsedMicroseconds.QuadPart = EndingTime.QuadPart - StartingTime.QuadPart;
            ElapsedMicroseconds.QuadPart *= 1000000;
            ElapsedMicroseconds.QuadPart /= Frequency.QuadPart;
            if (i == 0)
                time = (double)ElapsedMicroseconds.QuadPart / 1000;
            else {
                if (time > (double)ElapsedMicroseconds.QuadPart / 1000)
                    time = (double)ElapsedMicroseconds.QuadPart / 1000;
            }
        }
        sumtime += time;
    }
    printf("1000 loops,best of 3: %.3f ms per loop\n",sumtime/1000);

    free(Tm);
    free(tau);
    free(T);
}

Fortran f2py. Функцию можно использовать с Python.Версия от 2016 года (с функцией abs ()):

subroutine rec_fortran(tm,tau,alen,result)
    integer*8, intent(in) :: alen
    real*8, dimension(alen), intent(in) :: tm
    real*8, dimension(alen), intent(in) :: tau
    real*8, dimension(alen) :: res
    real*8, dimension(alen), intent(out) :: result

    res(1)=0
    do i=2,alen
        res(i) = tm(i) + (abs(res(i-1) - tm(i)))**(-tau(i))
    end do
    result=res    
end subroutine rec_fortran

Aronstef · Answer 4 · 18 мая 2017

Обновление: 21-10-2018 Я исправил свой ответ, основываясь на комментариях.

Возможно векторизовать операции над векторами, если вычисление не является рекурсивным.Поскольку рекурсивная операция зависит от предыдущего вычисленного значения, параллельная обработка операции невозможна.Следовательно, это не работает:

def calc_vect(Tm_, tau_):
    return Tm_[1:] - (Tm_[:-1] + Tm_[1:]) ** (-tau_[1:])

Поскольку (последовательная обработка / цикл) необходима, наилучшая производительность достигается при максимально возможном приближении к оптимизированному машинному коду, поэтому Numba и Cython - лучшие ответыздесь.

Подход Numba может быть реализован следующим образом:

init_string = """
from math import pow
import numpy as np
from numba import jit, float32

np.random.seed(0)
n = 100000
Tm = np.cumsum(np.random.uniform(0.1, 1, size=n).astype('float32'))
tau = np.random.uniform(-1, 0, size=n).astype('float32')

def calc_python(Tm_, tau_):
 tt = np.empty(len(Tm_))
 tt[0] = Tm_[0]
 for i in range(1, len(Tm_)):
     tt[i] = Tm_[i] - pow(tt[i-1] + Tm_[i], -tau_[i])
 return tt

@jit(float32[:](float32[:], float32[:]), nopython=False, nogil=True)
def calc_numba(Tm_, tau_):
  tt = np.empty(len(Tm_))
  tt[0] = Tm_[0]
  for i in range(1, len(Tm_)):
      tt[i] = Tm_[i] - pow(tt[i-1] + Tm_[i], -tau_[i])
  return tt
"""

import timeit
py_time = timeit.timeit('calc_python(Tm, tau)', init_string, number=100)
numba_time = timeit.timeit('calc_numba(Tm, tau)', init_string, number=100)
print("Python Solution: {}".format(py_time))
print("Numba Soltution: {}".format(numba_time))

Сравнение времени выполнения функций Python и Numba:

Python Solution: 54.58057559299999
Numba Soltution: 1.1389029540000024

parkus · Answer 5 · 13 ноября 2014

Чтобы основываться на ответе NPE, я согласен, что где-то должен быть цикл. Возможно, ваша цель - избежать накладных расходов, связанных с циклом Python for? В этом случае numpy.fromiter действительно выбивает цикл for, но только немного:

Используя очень простое рекурсивное соотношение,

x[i+1] = x[i] + 0.1

Я получаю

#FOR LOOP
def loopit(n):
     x = [0.0]
     for i in range(n-1): x.append(x[-1] + 0.1)
     return np.array(x)

#FROMITER
#define an iterator (a better way probably exists -- I'm a novice)
def it():
     x = 0.0
     while True:
         yield x
         x += 0.1

#use the iterator with np.fromiter
def fi_it(n):
     return np.fromiter(it(), np.float, n)

%timeit -n 100 loopit(100000)
#100 loops, best of 3: 31.7 ms per loop

%timeit -n 100 fi_it(100000)
#100 loops, best of 3: 18.6 ms per loop

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

def loopit(n):
     x = np.zeros(n)
     for i in range(n-1): x[i+1] = x[i] + 0.1
     return x

%timeit -n 100 loopit(100000)
#100 loops, best of 3: 50.1 ms per loop

Можно ли векторизовать рекурсивный расчет массива NumPy, где каждый элемент зависит от предыдущего?

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

Ответы [ 5 ]

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Можно ли векторизовать рекурсивный расчет массива NumPy, где каждый элемент зависит от предыдущего?

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

Ответы [ 5 ]

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Похожие темы