Нет скорости от Cython

Question

Я пытаюсь определить функцию, которая содержит внутренний цикл для моделирования интеграла.

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

Таким образом, я пытался заставить Cython работать, но я, должно быть, делаю серьезную ошибку (вероятно, многие из них!). Следуя документации Cython, я попытался напечатать мои переменные. После этого код работает так же медленно, как чистый Python. Это кажется странным.

Вот мой код:

import numpy as np 
cimport cython
cimport numpy as np
import minuit

data = np.genfromtxt('q6data.csv', usecols = np.arange(1, 24, 1), delimiter = ',')  

cdef int ns    = 1000                 # Number of simulation draws
cdef int K     = 5                    # Number of observed characteristics, including            constant
cdef int J     = len(data[:,1])       # Number of products, including outside
cdef double tol   = 0.0001            # Inner GMM loop tolerance
nu = np.random.normal(0, 1, (6, ns))  # ns random deviates

@cython.boundscheck(False)
@cython.wraparound(False)


def S(np.ndarray[double, ndim=1] delta, double s1, double s2, double s3, double s4,  double s5, double a):
    """Computes the simulated integrals, one for each good.
    Parameters: delta is an array of length J containing mean product specific utility levels
    Returns: Numpy array with length J."""
    cdef np.ndarray[double, ndim=2] mu_ij = np.dot((data[:,2:7]*np.array([s1, s2, s3, s4, s5])), nu[1:K+1,:])
    cdef np.ndarray[double, ndim=2] mu_y  = a * np.log(np.exp(data[:,21].reshape(J,1) +  data[:,22].reshape(J,1)*nu[0,:].reshape(1, ns)) - data[:,7].reshape(J,1))
    cdef np.ndarray[double, ndim=2] V = delta.reshape(J,1) + mu_ij + mu_y
    cdef np.ndarray[double, ndim=2] exp_vi = np.exp(V)
    cdef np.ndarray[double, ndim=2] P_i = (1.0 / np.sum(exp_vi[np.where(data[:,1] == 71)], 0)) * exp_vi[np.where(data[:,1] == 71)] 
    cdef int yrs = 19
    cdef int yr
    for yr in xrange(yrs):
        P_yr = (1.0 / np.sum(exp_vi[np.where(data[:,1]== (yr + 72))], 0)) *   exp_vi[np.where(data[:,1] == (yr + 72))]
        P_i  = np.concatenate((P_i, P_yr)) 
    cdef np.ndarray[double, ndim=1] S = np.zeros(dtype = "d", shape = J)
    cdef int j
    for j in xrange(ns):
        S += P_i[:,j]
    return (1.0 / ns) * S

def d_infty(np.ndarray[double, ndim=1] x, np.ndarray[double, ndim=1] y):
    """Sup norm."""
    return np.max(np.abs(x - y)) 

def T(np.ndarray[double, ndim=1] delta_exp, double s1, double s2, double s3, double s4,  double s5, double a):
    """The contraction operator.  This function takes the parameters and the exponential
    of the starting value of delta and returns the fixed point.""" 
    cdef int iter = 0
    cdef int maxiter = 200
    cdef int i
    for i in xrange(maxiter): 
        delta1_exp = delta_exp * (data[:, 8] / S(np.log(delta_exp), s1, s2, s3, s4, s5, a))                    
        print i
        if d_infty(delta_exp, delta1_exp) < tol:                                       
            break
        delta_exp = delta1_exp
    return np.log(delta1_exp)


def Q(double s1, double s2, double s3, double s4, double s5, double a):
    """GMM objective function."""  
    cdef np.ndarray[double, ndim=1] delta0_exp = np.exp(data[:,10])                                                     
    cdef np.ndarray[double, ndim=1] delta1 = T(delta0_exp, s1, s2, s3, s4, s5, a)
    delta1[np.where(data[:,10]==0)] = np.zeros(len(np.where(data[:,10]==0)))            
    cdef np.ndarray[double, ndim=1] xi =  delta1 - (np.dot(data[:,2:7],   np.linalg.lstsq(data[:,2:7], delta1)[0]))   
    cdef np.ndarray[double, ndim=2] g_J = xi.reshape(J,1) * data[:,11:21]
    cdef np.ndarray[double, ndim=1] G_J = (1.0 / J) * np.sum(g_J, 0) 
    return np.sqrt(np.dot(G_J, G_J))

Я профилировал код, и, похоже, это функция S, встроенный симулятор, который убивает производительность. Во всяком случае, я бы ожидал хоть какой-то прирост скорости от набора переменных. Поскольку это не принесло никакой пользы, я склонен полагать, что я делаю некоторые фундаментальные ошибки.

Кто-нибудь видит явную ошибку в коде Cython, которая может привести к такому результату?

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

Answer 1

Cython может создать html-файл, чтобы помочь с этим:

cython -a MODULE.py

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

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

Без этого легко пропустить некоторые статические объявления типов переменных, вызовов функций или операторов.(например, оператор индексации x [y] по-прежнему является полностью динамической операцией Python, если вы не объявите иначе)

Answer 2

Cython не предлагает автоматическое повышение производительности, вы должны знать его внутренние компоненты и проверять сгенерированный код C.

В частности, если вы хотите улучшить производительность циклов, вам следует избегать вызова в них функций Python, которые вы часто делаете в этом случае (все вызовы np. - это вызовы Python, срезы и, возможно, другие вещи). ).

См. на этой странице общие рекомендации по оптимизации производительности с помощью Cython (ключ -a действительно удобен при оптимизации) и этот для подробностей при оптимизации кода Numpy.

Answer 3

Вы можете определенно ускорить свой код, используя больше возможностей Numpy.

Например:

cdef np.ndarray[double, ndim=1] S = np.zeros(dtype = "d", shape = J)
cdef int j
for j in xrange(ns):
    S += P_i[:,j]

будет намного быстрее и удобочитаемее, как

S = P_i.sum(axis=1)

Вы также повторяете некоторые вычисления, которые, таким образом, занимают вдвое больше времени, чем необходимо.Например,

np.where(data[:,1]==(yr + 72))

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

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

Answer 4

Профилировщик поможет вам выяснить, какая часть медленная?Мне нравится запускать программы с использованием стандартного библиотечного профилировщика:

python -O -m cProfile -o profile.out MYAPP.py

, а затем просматривать вывод этого в графическом интерфейсе RunSnakeRun:

runsnake profile.out

Отсюда можно установить RunSnakeRun: http://www.vrplumber.com/programming/runsnakerun/

Answer 5

Разделение данных только , после вашего комментария 28 ноября:

import sys
import time
import numpy as np

def splitdata( data, n, start=1971 ):
    """ split data into n pieces, where col 1 == start .. start + n """
        # not fancy, fast enough for small n
    split = n * [None]
    for j in range(n):
        split[j] = data[ data[:,1] == start + j ]
    return split  # [ arrays: col1 0, col1 1 ... ]

#...........................................................................
N = 2237
ncol = 21
start = 1971
n = 20
seed = 1
exec "\n".join( sys.argv[1:] )  # run this.py N= ...
np.set_printoptions( 2, threshold=100, suppress=True )  # .2f
np.random.seed(seed)

print "N=%d  ncol=%d  n=%d" % (N, ncol, n)
data = np.random.uniform( start, start + n, (N,ncol) )
data[:,1] = data[:,1].round()
t0 = time.time()

split = splitdata( data, n, start )  # n pieces

print "time: %.0f us splitdata" % ((time.time() - t0) * 1e6)
for y, yeardata in enumerate(split):
    print "%d  %d  %g" % (start + y, len(yeardata), yeardata[:,0].sum())

->

time: 27632 us splitdata  # old mac ppc
1971  69  136638
1972  138  273292
...

Answer 6

Следуя приведенным здесь советам, я потратил больше времени на профилирование приведенного выше кода. Надеюсь, я кое-что почистил ...

Я немного профилировал код и лучше представляю, какие части самые медленные. Я дополнительно определил

X = data[:, 2:7]
m_y = data[:, 21].reshape(J,1)
sigma_y = 1.0
price = data[:, 7].reshape(J, 1)
shares_data = data[:,8]

Тогда это следующие строки, которые съедают большую часть общего времени.

mu_ij = np.dot((X*np.array([s1, s2, s3, s4, s5])), nu[1:K+1,:])
mu_y  = a * np.log(np.exp(m_y + sigma_y*nu[0,:].reshape(1,ns)) - price)
V = delta.reshape(J,1) + mu_ij + mu_y
exp_vi = np.exp(V)
P_i = (1.0 / np.sum(exp_vi[np.where(data[:,1]==71)], 0)) *  exp_vi[np.where(data[:,1]==71)] 
for yr in xarange(19):
    P_yr = (1.0 / np.sum(exp_vi[np.where(data[:,1]==yr)], 0)) * exp_vi[np.where(data[:,1]==yr)]
P_i  = np.concatenate((P_i, P_yr))

У меня сложилось впечатление, что это слишком громоздкий способ достижения моей цели. Я надеялся, что кто-нибудь сможет дать совет о том, как ускорить работу этих линий. Может быть, есть какие-то возможности Numpy, которые мне не хватает? Если эта проблема недостаточно четко определена, чтобы вы были полезны, я был бы рад предоставить более подробную информацию о контексте моей проблемы. Спасибо!

Answer 7

«Фундаментальная ошибка» заключается в том, что вы ожидаете хорошей производительности в длинных циклах от Python.Это интерпретируемый язык, и переключение между реализациями и ctyping ничего не делает для этого.Есть несколько числовых библиотек Python для быстрых вычислений, написанных в основном на C. Например, если вы уже используете numpy для массивов, почему бы вам не пойти дальше и использовать scipy для продвинутой математики?Это увеличит скорость чтения и .