Блог

Главная - Вопросы по программированию - Сравнение scipy.stats.t.sf и GSL с использованием Cython

Сравнение scipy.stats.t.sf и GSL с использованием Cython

#python #numpy #scipy #cython #gsl

#питон #numpy #сципи #cython #gsl #python #scipy

Вопрос:

Я хотел бы вычислить значения p большого 2D-массива значений numpy t. Однако это занимает много времени, и я хотел бы улучшить его скорость. Я попробовал использовать GSL. Хотя один gsl_cdf_tdist_P намного быстрее, чем scipy.stats.t.sf, при итерации по ndarray процесс выполняется очень медленно. Я хотел бы получить помощь, чтобы улучшить это. Смотрите код ниже.

GSL_Test.pyx

 import cython
cimport cython

import numpy
cimport numpy
from cython_gsl cimport *

DTYPE = numpy.float32
ctypedef numpy.float32_t DTYPE_t

cdef get_gsl_p(double t, double nu):
    return (1 - gsl_cdf_tdist_P(t, nu)) * 2

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cdef get_gsl_p_for_2D_matrix(numpy.ndarray[DTYPE_t, ndim=2] t_matrix, int n):

    cdef unsigned int rows = t_matrix.shape[0]
    cdef numpy.ndarray[DTYPE_t, ndim=2] out = numpy.zeros((rows, rows), dtype='float32')
    cdef unsigned int row, col

    for row in range(rows):
        for col in range(rows):
            out[row, col] = get_gsl_p(t_matrix[row, col], n-2)
    return out

def get_gsl_p_for_2D_matrix_def(numpy.ndarray[DTYPE_t, ndim=2] t_matrix, int n):
    return get_gsl_p_for_2D_matrix(t_matrix, n)

ipython

 import GSL_Test
import numpy
import scipy.stats
a = numpy.random.rand(3544, 3544).astype('float32')
%timeit -n 1 GSL_Test.get_gsl_p_for_2D_matrix(a, 25)
1 loop, best of 3: 7.87 s per loop
%timeit -n 1 scipy.stats.t.sf(a, 25)*2
1 loop, best of 3: 4.66 s per loop

ОБНОВЛЕНИЕ: добавив объявления cdef, я смог сократить время вычислений, но по-прежнему не ниже, чем scipy. Я изменил код, чтобы иметь объявления cdef. 

 %timeit -n 1 GSL_Test.get_gsl_p_for_2D_matrix_def(a, 25)
1 loop, best of 3: 6.73 s per loop

Комментарии:

1. О чем вы знаете stats.t.sf ? Как вы думаете, почему cython версия может быть быстрее? Я ищу такие вещи, как итерации, на которые можно перевести c .

2. scipy.stats.t.sf и другие методы имеют некоторые накладные расходы python, в основном для проверки ввода. Если вам нужна скорость, вы можете напрямую использовать функции scipy.special, которые написаны либо на Fortran, либо на C. Я сомневаюсь, что можно превзойти scipy.special на сколько-нибудь значительную величину, если не существует более эффективного алгоритма или алгоритма, который обменивает скорость на точность.

3. Нет никакого преимущества в том, чтобы иметь отдельную def cdef функцию and для get_gsl_p_for_2D_matrix . Вы можете получить немного информации, указав возвращаемый тип для get_gsl_p .

Ответ №1:

Вы можете получить небольшой выигрыш в производительности raw, используя специальную функцию raw вместо stats.t.sf . Глядя на источник, вы обнаружите (https://github.com/scipy/scipy/blob/master/scipy/stats/_continuous_distns.py#L3849 )

 def _sf(self, x, df):
      return sc.stdtr(df, -x)

Чтобы вы могли использовать stdtr напрямую:

 np.random.seed(1234)
x = np.random.random((3740, 374))
t1 = stats.t.sf(x, 25)
t2 = stdtr(25, -x)

1 loop, best of 3: 653 ms per loop
1 loop, best of 3: 562 ms per loop

Если вы обращаетесь к cython, типизированный синтаксис memoryview часто дает вам более быстрый код, чем старый синтаксис ndarray:

 from scipy.special.cython_special cimport stdtr
from numpy cimport npy_intp
import numpy as np

def tsf(double [:, ::1] x, int df=25):
    cdef double[:, ::1] out = np.empty_like(x)
    cdef npy_intp i, j
    cdef double tmp, xx

    for i in range(x.shape[0]):
        for j in range(x.shape[1]):
            xx = x[i, j]
            out[i, j] = stdtr(df, -xx)

    return np.asarray(out)

Здесь я также использую cython_special интерфейс, который доступен только в версии scipy для разработчиков (http://scipy.github.io/devdocs/special.cython_special.html#module-scipy.special.cython_special ), но вы можете использовать GSL, если хотите.

Наконец, если вы подозреваете узкое место в итерациях, не забудьте проверить вывод cython -a , чтобы увидеть, есть ли какие-либо накладные расходы python в горячих циклах.