Оптимизация моего кода Cython/Numpy? Прирост производительности пока только 30%

Есть ли что-то, что я забыл сделать здесь, чтобы немного ускорить процесс? Я пытаюсь реализовать алгоритм, описанный в книге Tuning Timbre Spectrum Scale. Кроме того, если ничего не помогает, есть ли способ просто написать эту часть кода на C, а затем вызвать ее из python?

import numpy as np
cimport numpy as np

# DTYPE = np.float
ctypedef np.float_t DTYPE_t

np.seterr(divide='raise', over='raise', under='ignore', invalid='raise')

"""
I define a timbre as the following 2d numpy array:
[[f0, a0], [f1, a1], [f2, a2]...] where f describes the frequency
of the given partial and a is its amplitude from 0 to 1. Phase is ignored.
"""

#Test Timbre
# cdef np.ndarray[DTYPE_t,ndim=2] t1 = np.array( [[440,1],[880,.5],[(440*3),.333]])

# Calculates the inherent dissonance of one timbres of the above form
# using the diss2Partials function
cdef DTYPE_t diss1Timbre(np.ndarray[DTYPE_t,ndim=2] t):
    cdef DTYPE_t runningDiss1
    runningDiss1 = 0.0
    cdef unsigned int len = np.shape(t)[0]
    cdef unsigned int i
    cdef unsigned int j
    for i from 0 <= i < len:
        for j from i+1 <= j < len:
            runningDiss1 += diss2Partials(t[i], t[j])
    return runningDiss1

# Calculates the dissonance between two timbres of the above form 
cdef DTYPE_t diss2Timbres(np.ndarray[DTYPE_t,ndim=2] t1, np.ndarray[DTYPE_t,ndim=2] t2):
    cdef DTYPE_t runningDiss2
    runningDiss2 = 0.0
    cdef unsigned int len1 = np.shape(t1)[0]
    cdef unsigned int len2 = np.shape(t2)[0]
    runningDiss2 += diss1Timbre(t1)
    runningDiss2 += diss1Timbre(t2)
    cdef unsigned int i1
    cdef unsigned int i2
    for i1 from 0 <= i1 < len1:
        for i2 from 0 <= i2 < len2:
            runningDiss2 += diss2Partials(t1[i1], t2[i2])
    return runningDiss2

cdef inline DTYPE_t float_min(DTYPE_t a, DTYPE_t b): return a if a <= b else b

# Calculates the dissonance of two partials of the form [f,a]
cdef DTYPE_t diss2Partials(np.ndarray[DTYPE_t,ndim=1] p1, np.ndarray[DTYPE_t,ndim=1] p2):
    cdef DTYPE_t f1 = p1[0]
    cdef DTYPE_t f2 = p2[0]
    cdef DTYPE_t a1 = abs(p1[1])
    cdef DTYPE_t a2 = abs(p2[1])

    # In order to insure that f2 > f1:
    if (f2 < f1):
        (f1,f2,a1,a2) = (f2,f1,a2,a1)

    # Constants of the dissonance curves
    cdef DTYPE_t _xStar
    _xStar = 0.24
    cdef DTYPE_t _s1
    _s1 = 0.021
    cdef DTYPE_t _s2
    _s2 = 19
    cdef DTYPE_t _b1
    _b1 = 3.5
    cdef DTYPE_t _b2
    _b2 = 5.75

    cdef DTYPE_t a = float_min(a1,a2)
    cdef DTYPE_t s = _xStar/(_s1*f1 + _s2)
    return (a * (np.exp(-_b1*s*(f2-f1)) - np.exp(-_b2*s*(f2-f1)) ) )

cpdef dissTimbreScale(np.ndarray[DTYPE_t,ndim=2] t,np.ndarray[DTYPE_t,ndim=1] s):
    cdef DTYPE_t currDiss
    currDiss = 0.0;
    cdef unsigned int i
    for i from 0 <= i < s.size:
        currDiss += diss2Timbres(t, transpose(t,s[i]))
    return currDiss

cdef np.ndarray[DTYPE_t,ndim=2] transpose(np.ndarray[DTYPE_t,ndim=2] t, DTYPE_t ratio):
    return np.dot(t, np.array([[ratio,0],[0,1]]))

Ссылка на код: Cython Code


person Chironex    schedule 17.03.2011    source источник
comment
например: docs.cython.org/src/tutorial/external.html и docs.cython.org/src/tutorial/clibraries.html?   -  person reve_etrange    schedule 17.03.2011
comment
Нет, я использую Cython только для ускорения алгоритма, изначально написанного на чистом Python.   -  person Chironex    schedule 17.03.2011
comment
Для python -> C см. SO с-китоном   -  person denis    schedule 17.03.2011


Ответы (3)


Вот некоторые вещи, которые я заметил:

  1. Используйте t1.shape[0] вместо np.shape(t1)[0] и так далее в других местах.
  2. Не используйте len в качестве переменной, потому что это встроенная функция в Python (не для скорости, а для хорошей практики). Используйте L или что-то в этом роде.
  3. Не передавайте двухэлементные массивы в функции, если вам это действительно не нужно. Cython проверяет буфер каждый раз, когда вы передаете массив. Итак, при использовании diss2Partials(t[i], t[j]) вместо этого сделайте diss2Partials(t[i,0], t[i,1], t[j,0], t[j,1]) и соответствующим образом переопределите diss2Partials.
  4. Не используйте abs или, по крайней мере, не Python. Необходимо преобразовать ваш C-двойник в число с плавающей запятой Python, вызвать функцию abs, а затем преобразовать обратно в C-двойник. Вероятно, было бы лучше сделать встроенную функцию, как вы сделали с float_min.
  5. Вызов np.exp аналогичен использованию abs. Измените np.exp на exp и добавьте from libc.math cimport exp к вашему импорту вверху.
  6. Полностью избавьтесь от функции transpose. np.dot действительно замедляет работу, но в любом случае здесь нет необходимости в умножении матриц. Перепишите свою функцию dissTimbreScale, чтобы создать пустую матрицу, скажем, t2. Перед текущим циклом установите второй столбец t2 равным второму столбцу t (предпочтительно с использованием цикла, но вы, вероятно, можете обойтись здесь операцией Numpy). Затем внутри текущего цикла добавьте цикл, который устанавливает первый столбец t2 равным первому столбцу t, умноженному на s[i]. Вот что на самом деле делало ваше матричное умножение. Затем просто передайте t2 в качестве второго параметра функции diss2Timbres вместо возвращаемого функцией transpose.

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

person Justin Peel    schedule 24.03.2011
comment
Я заменил Numpy exps на Python exps, но вау! Использование C exps приводит к хорошему 2-кратному ускорению! Теперь рассмотрим №3 и №6... - person Chironex; 25.03.2011
comment
Хорошо, реализовали #3 и #6 и поэкспериментировали с тем, как код обрабатывает несортированные списки (вместо того, чтобы сортировать и затем передавать их, он оставляет их несортированными и ожидает, что следующая функция будет сравнивать каждую пару элементов, чтобы использовать их в правильный порядок)... в любом случае, большое спасибо за комментарий, он сделал мой код в 2,5 раза быстрее (после того, как я уже сделал его в 15 раз быстрее)! - person Chironex; 25.03.2011
comment
Хммм... может мне стоит добавить функцию быстрой сортировки, чтобы мне не приходилось делать все эти сравнения... мда - person Chironex; 25.03.2011

В вашем коде:

for i from 0 <= i < len:
    for j from i+1 <= j < len:
        runningDiss1 += diss2Partials(t[i], t[j])
return runningDiss1

проверка границ выполняется для каждого поиска массива, используйте декоратор @cython.boundscheck(False) перед функцией, а затем приводите к типу unsigned int перед использованием i и j в качестве индексов. Дополнительную информацию см. в руководстве по Cython для Numpy.

person highBandWidth    schedule 17.03.2011

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

Я сравнил Python/Cython и Numexpr для одной из моих функций (ссылка на SO). В зависимости от размера массива numexpr превосходил как Cython, так и Fortran.

ПРИМЕЧАНИЕ. Только что понял, что этот пост действительно старый...

person Moritz    schedule 13.02.2015