Пользователь
0,0
рейтинг
30 сентября 2014 в 17:21

Разработка → И еще раз о GIL в Python из песочницы

Предисловие


Область, в которой мне повезло работать, называется вычислительная электрофизиология сердца. Физиология сердечной деятельности определяется электрическими процессами, происходящими на уровне отдельных клеток миокарда. Эти электрические процессы создают электрическое поле, которое достаточно легко измерить. Более того оно очень неплохо описывается в рамках математических моделей электростатики. Тут и возникает уникальная возможность строго математически описать работу сердца, а значит — и усовершенствовать методы лечения многих сердечных заболеваний.

За время работы в этой области у меня накопился некоторый опыт использования различных вычислительных технологий. На некоторые вопросы, которые могут быть интересны не только мне, я постараюсь отвечать в рамках этой публикации.

Кратко о Scientific Python


Начиная еще с первых курсов университета, я пытался найти идеальный инструмент для быстрой разработки численных алгоритмов. Если отбросить ряд откровенно маргинальных технологий, я курсировал между C++ и MATLAB. Это продолжалось до тех пор, пока я не открыл для себя Scientific Python [1].

Scientific Python представляет собой набор библиотек языка Python для научных вычислений и научной визуализации. В своей работе я использую следующие пакеты, которые покрывают примерно 90% моих потребностей:
Название Описание
NumPy Одна из базовых библиотек, позволяет работать с многомерными массивами как с едиными объектами в MATLAB стиле.
Включает реализацию основных процедур линейной алгебры, преобразование Фурье, работу со случайными числами и др.
SciPy Расширение NumPy, включает реализацию методов оптимизации, работу с разряженными матрицами, статистику и др.
Pandas Отдельный пакет для анализа многомерных данных и статистики.
SymPy Пакет символьной математики.
Matplotlib Двумерная графика.
Mayavi2 Трехмерная графика на основе VTK.
Spyder Удобная IDE для интерактивной разработки математических алгоритмов.
В Scientific Python я нашел для себя великолепный баланс между удобной высокоуровневой абстракцией для быстрой разработки численных алгоритмов и современным развитым языком. Но, как известно, не бывает идеальных инструментов. И одна из достаточно критических проблем в Python — это проблема параллельных вычислений.

Проблемы параллельных вычислений в Python.


Под параллельными вычислениями в этой статье я буду понимать SMP — симметричный мультипроцессинг с общей памятью. Вопросов использования CUDA и систем с раздельной памятью (чаще всего используется стандарт MPI) касаться не буду.

Проблема заключается в GIL. GIL (Global Interpreter Lock) — это блокировка (mutex), которая не позволяет нескольким потокам выполнить один и тот же байткод. Эта блокировка, к сожалению, является необходимой, так как система управления памятью в CPython не является потокобезопасной. Да, GIL это не проблема языка Python, а проблема реализации интерпретатора CPython. Но, к сожалению, остальные реализации Python не слишком приспособлены для создания быстрых численных алгоритмов.

К счастью, в настоящее время существует несколько способов решения проблем GIL. Рассмотрим их.

Тестовая задача


Даны два набора по N векторов: P={p1,p2,…,pN} и Q={q1,q2,…,qN} в трехмерном евклидовом пространстве. Необходимо построить матрицу R размерностью N x N, каждый элемент ri,j которой вычисляется по формуле:


Грубо говоря, нужно вычислить матрицу, использующую попарные расстояния между всеми векторами. Эта матрица достаточно часто используется в реальных расчетах, например, при RBF интерполяции или решении дифуров в чп методом интегральных уравнений.

В тестовых экспериментах количество векторов N = 5000. Для вычислений использовался процессор с 4 ядрами. Результаты получены по среднему времени из 10 запусков.

Полную реализацию тестовых задач можно поглядеть на GitHub [2].
Правильное замечание в комментариях от "@chersaya". Данная тестовая задача используется здесь в качестве примера. Если нужно действительно вычислить попарные расстояния, правильнее использовать функцию scipy.spatial.distance.cdist.

Параллельная реализация на C++


Для сравнения эффективности параллельных вычислений на Python, я реализовал эту задачу на C++. Код основной функции выглядит следующий образом.

Однопроцессорная реализация:

//! Single thread matrix R calculation
void spGetR(vector<Vector3D> & p, vector<Vector3D> & q, MatrixMN & R)
{
    for (int i = 0; i < p.size(); i++) {
        Vector3D & a = p[i];
        for (int j = 0; j < q.size(); j++) {
            Vector3D & b = q[j];

            Vector3D r = b - a;
            R(i, j) = 1 / (1 + sqrt(r * r));
        }
    }
}

Многопроцессорная реализация:

//! OpenMP matrix R calculations
void mpGetR(vector<Vector3D> & p, vector<Vector3D> & q, MatrixMN & R)
{
    #pragma omp parallel for
    for (int i = 0; i < p.size(); i++) {
        Vector3D & a = p[i];
        for (int j = 0; j < q.size(); j++) {
            Vector3D & b = q[j];

            Vector3D r = b - a;
            R(i, j) = 1 / (1 + sqrt(r * r));
        }
    }
}

Что здесь интересного? Ну прежде всего я использовал отдельный класс Vector3D для представления вектора в трехмерном пространстве. Перегруженный оператор "*" в этом классе имеет смысл скалярного произведения. Для представления набора векторов я использовал std::vector. Для параллельных вычислений использовалась технология OpenMP. Для параллелизации алгоритма достаточно использовать директиву "#pragma omp parallel for".

Результаты:
Однопроцессорный С++ 224 ms
Многопроцессорный C++ 65 ms
Ускорение в 3.45 раза при параллельном расчете я считаю вполне неплохим для четырехядерного процессора.

Параллельные реализации на Python


1.Наивная реализация на чистом Python

В этом тесте хотелось проверить сколько будет решаться задача на чистом Python без использования каких-либо специальных пакетов.

Код решения:

def sppyGetR(p, q):
    R = np.empty((p.shape[0], q.shape[1]))    
    
    nP = p.shape[0]
    nQ = q.shape[1]
    for i in xrange(nP):
        for j in xrange(nQ):
            rx = p[i, 0] - q[0, j]
            ry = p[i, 1] - q[1, j]
            rz = p[i, 2] - q[2, j]

            R[i, j] = 1 / (1 + sqrt(rx * rx + ry * ry + rz * rz))
    return R

Здесь p, q – входные данные в формате NumPy массивов размерностями (N, 3) и (3, N). А дальше идет честный цикл на Python, вычисляющий элементы матрицы R.

Результаты:
Однопроцессорный Python 57 386 ms
Да, да, именно 57 тысяч миллисекунд. Где-то в 256 раз медленнее однопроцессорного C++. В общем, это совсем не вариант для численных расчетов.

2 Однопроцессорный NumPy

Вообще, для вычислений на Python с использованием NumPy иногда можно вообще не задумываться о параллельности. Так, например, процедура умножения двух матриц на NumPy будет в итоге все-равно выполняться с использованием низкоуровневых высокоэффективных библиотек линейной алгебры на C++ (MKL или ATLAS). Но, к сожалению, это верно лишь для наиболее типовых операций и не работает в общем случае. Наша тестовая задача, к сожалению, будет выполняться последовательно.

Код решения следующий:

def spnpGetR(p, q):    
    Rx = p[:, 0:1] - q[0:1]
    Ry = p[:, 1:2] - q[1:2]
    Rz = p[:, 2:3] - q[2:3]
         
    R = 1 / (1 + np.sqrt(Rx * Rx + Ry * Ry + Rz * Rz))
       
    return R


Всего 4 строчки и никаких циклов! Вот за это я и люблю NumPy.

Результаты:
Однопроцессорный NumPy 973 ms
Примерно в 4.3 раза медленнее однопроцессорного C++. Вот это уже совсем неплохой результат. Для подавляющего большинства расчетов этой производительности вполне хватает. Но это все пока однопроцессорные результаты. Идем дальше к мультипроцессингу.

3 Многопроцессорный NumPy

В качестве решения проблем с GIL традиционно предлагается использовать несколько независимых процессов выполнения вместо нескольких потоков выполнения. Все бы хорошо, но есть проблема. Каждый процесс обладает независимой памятью, и нам необходимо в каждый процесс передавать матрицу результатов. Для решения этой проблемы в Python multiprocessing вводится класс RawArray, предоставляющий возможность разделить один массив данных между процессами. Не знаю точно, что лежит в основе RawArray. Мне кажется, что это memory mapped files.

Код решения следующий:

def mpnpGetR_worker(job):
    start, stop = job
    p = np.reshape(np.frombuffer(mp_share.p), (-1, 3))
    q = np.reshape(np.frombuffer(mp_share.q), (3, -1))
    R = np.reshape(np.frombuffer(mp_share.R), (p.shape[0], q.shape[1]))

    Rx = p[start:stop, 0:1] - q[0:1]
    Ry = p[start:stop, 1:2] - q[1:2]
    Rz = p[start:stop, 2:3] - q[2:3]

    R[start:stop, :] = 1 / (1 + np.sqrt(Rx * Rx + Ry * Ry + Rz * Rz))
        
def mpnpGetR(p, q):
    nP, nQ = p.shape[0], q.shape[1]    
    
    sh_p = mp.RawArray(ctypes.c_double, p.ravel())
    sh_q = mp.RawArray(ctypes.c_double, q.ravel())
    sh_R = mp.RawArray(ctypes.c_double, nP * nQ)
    
    nCPU = 4
    jobs = utils.generateJobs(nP, nCPU)
    
    pool = mp.Pool(processes=nCPU, initializer=mp_init, initargs=(sh_p, sh_q, sh_R))
    pool.map(mpnpGetR_worker, jobs, chunksize=1)

    R = np.reshape(np.frombuffer(sh_R), (nP, nQ))
    return R

Мы создаем разделенные массивы для входных данных и выходной матрицы, создаем пул процессов по числу ядер, разбиваем задачу на подзадачи и решаем параллельно.

Результаты:
Многопроцессорный NumPy 795 ms
Да, быстрее однопроцессорного варианта, но всего в 1.22 раза. С ростом числа N эффективность решения растет. Но, в целом и общем, наша тестовая задача не слишком приспособлена для решения в рамках множества независимых процессов с независимой памятью. Хотя для других задач такой вариант может быть вполне эффективным.

На этом известные мне решения для параллельного программирования с использование только Python закончились. Далее, как бы нам не хотелось, для освобождения от GIL придется спускаться на уровень C++. Но этот не так страшно, как кажется.

4 Cython

Cython [3] — это расширение языка Python, позволяющее внедрять инструкции на языке C в код на Python. Таким образом, мы можем взять код на Python и добавлением нескольких инструкций значительно ускорить узкие в плане производительности места. Cython модули преобразуются в код на C и далее компилируются в Python модули. Код решения нашей задачи на Cython следующий:

Однопроцессорный Cython:

@cython.boundscheck(False)
@cython.wraparound(False)
def spcyGetR(pp, pq):
    pR = np.empty((pp.shape[0], pq.shape[1]))    
    
    cdef int i, j, k
    cdef int nP = pp.shape[0]
    cdef int nQ = pq.shape[1]
    cdef double[:, :] p = pp
    cdef double[:, :] q = pq
    cdef double[:, :] R = pR
    cdef double rx, ry, rz
    
    with nogil:
        for i in xrange(nP):
            for j in xrange(nQ):
                rx = p[i, 0] - q[0, j]
                ry = p[i, 1] - q[1, j]
                rz = p[i, 2] - q[2, j]

                R[i, j] = 1 / (1 + sqrt(rx * rx + ry * ry + rz * rz))
    return R

Многопроцессорный Cython:
@cython.boundscheck(False)
@cython.wraparound(False)
def mpcyGetR(pp, pq):
    pR = np.empty((pp.shape[0], pq.shape[1]))    
    
    cdef int i, j, k
    cdef int nP = pp.shape[0]
    cdef int nQ = pq.shape[1]
    cdef double[:, :] p = pp
    cdef double[:, :] q = pq
    cdef double[:, :] R = pR
    cdef double rx, ry, rz
    
    with nogil, parallel():
        for i in prange(nP, schedule='guided'):
            for j in xrange(nQ):
                rx = p[i, 0] - q[0, j]
                ry = p[i, 1] - q[1, j]
                rz = p[i, 2] - q[2, j]
                
                R[i, j] = 1 / (1 + sqrt(rx * rx + ry * ry + rz * rz))
                
    return R

Если сравнить данный код с реализацией на чистом Python, то все, что нам пришлось сделать, это всего лишь указать типы для используемых переменных. GIL отпускается одной строчкой. Параллельный цикл организуется всего лишь инструкцией prange вместо xrange. На мой взгляд, вполне несложно и красиво!

Результаты:
Однопроцессорный Cython 255 ms
Многопроцессорный Cython 75 ms
Вау! Время исполнения почти совпадает с временем исполнения на C++. Отставание примерно в 1.1 раз как в однопроцессорном, так и в многопроцессорном вариантах практически незаметно на реальных задачах.

5 Numba

Numba [4] достаточно новая библиотека, находится в активном развитии. Идея здесь примерно такая же, что и в Cython — попытка спуститься на уровень C++ в коде на Python. Но идея реализована существенно элегантнее.

Numba основана на LLVM компиляторах, которые позволяют производить компиляцию непосредственно в процессе исполнения программы (JIT компиляция). Например, для компилирования любой процедуры на Python достаточно всего лишь добавить аннотацию «jit». Более того аннотации позволяют указывать типы входных/выходных данных, что делает JIT-компиляцию существенно более эффективной.
Код реализации задачи следующий.

Однопроцессорный Numba:

@jit(double[:, :](double[:, :], double[:, :]))
def spnbGetR(p, q):
    nP = p.shape[0]
    nQ = q.shape[1]
    
    R = np.empty((nP, nQ))    
    
    for i in xrange(nP):
        for j in xrange(nQ):
            rx = p[i, 0] - q[0, j]
            ry = p[i, 1] - q[1, j]
            rz = p[i, 2] - q[2, j]

            R[i, j] = 1 / (1 + sqrt(rx * rx + ry * ry + rz * rz))
    return R

Многопроцессорный Numba:
def makeWorker():
    savethread = pythonapi.PyEval_SaveThread
    savethread.argtypes = []
    savethread.restype = c_void_p

    restorethread = pythonapi.PyEval_RestoreThread
    restorethread.argtypes = [c_void_p]
    restorethread.restype = None

    def worker(p, q, R, job):
        threadstate = savethread()
    
        nQ = q.shape[1]
        for i in xrange(job[0], job[1]):
            for j in xrange(nQ):
                rx = p[i, 0] - q[0, j]
                ry = p[i, 1] - q[1, j]
                rz = p[i, 2] - q[2, j]

                R[i, j] = 1 / (1 + sqrt(rx * rx + ry * ry + rz * rz))
    
        restorethread(threadstate)
        
    signature = void(double[:, :], double[:, :], double[:, :], int64[:])
    worker_ext = jit(signature, nopython=True)(worker)
    
    return worker_ext

def mpnbGetR(p, q):
    nP, nQ = p.shape[0], q.shape[1]
    R = np.empty((nP, nQ))

    nCPU = utils.getCPUCount()
    jobs = utils.generateJobs(nP, nCPU) 
    worker_ext = makeWorker()

    threads = [threading.Thread(target=worker_ext, args=(p, q, R, job)) for job in jobs]
    for thread in threads:
        thread.start()
    for thread in threads:
        thread.join()

    return R

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

Результаты выполнения:
Однопроцессорный Numba 359 ms
Многопроцессорный Numba 180 ms
Слегка хуже, чем Cython, но результаты все-равно очень достойные! А само решение крайне элегантное.

Выводы


Результаты я хочу проиллюстрировать следующими диаграммами:

image

Рис. 1. Результаты однопроцессорных вычислений

image

Рис. 2. Результаты многопроцессорных вычислений

Мне кажется, что проблемы GIL в Python для численных расчетов практически преодолены. Пока в качестве технологии параллельных вычислений я бы рекомендовал Cython. Но очень бы внимательно пригляделся бы к Numba.

Ссылки


[1] Scientific Python: scipy.org
[2] Полные исходные коды тестов: github.com/alec-kalinin/open-nuance
[3] Cython: cython.org
[4] Numba: numba.pydata.org

P.S. В комментариях "@chersaya" правильно указал еще один способ параллельных вычислений. Это использование библиотеки numexpr. Numexpr использует собственную виртуальную машину, написанную на C, и собственный JIT компилятор. Это позволяет ему принимать простые математические выражения в виде строки, компилировать и быстро вычислять.

Пример использования:
import numpy as np
import numexpr as ne

a = np.arange(1e6)
b = np.arange(1e6)

result = ne.evaluate("sin(a) + arcsinh(a/b)")
Alexander Kalinin @alec_kalinin
карма
30,0
рейтинг 0,0
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Спецпроект

Самое читаемое Разработка

Комментарии (47)

  • –9
    И вот почему в Python 3 вместо решения одной проблемы — GIL — создали дополнительные в виде юникода, бинарей и др?
    • +10
      И вот почему некоторые не видят разницы между языком и его интерпретаторами?
      • +1
        И почему некоторые считают, что язык, у которого только одна реализация, и на особенности которой уже заточена вся инфраструктура, не равен его интерпретатору?

        Когда-нибудь обязательно настанет счастье, но сейчас…
        • +3
          Основные реализации — CPython, PyPy, IronPython, Jython
          • 0
            Да, сейчас стало чуть легче, и IronPython более-менее исправляет ситуацию. Но неслабый пласт модулей прибиты к CPython гвоздями.
            Тот же numpy, matplotlib, pygtk…
      • +29
        Люди, которые не понимают, зачем в Py3 unicode стал дефолтным строковым типом, раньше говнокодили страницы на cp1251 в пхп.
        • +4
          Не только «говнокодили страницы на cp1251 в пхп», но ещё писали под DOS, например (когда Python собственно и не было).

          Реально Py3 и Py2 — это разные языки. И продолжение Py2 это скорее Pyston, а не Py3. А кому нужен дефолтный юникод, скорее будет уже выбирать между Py3 и Go.
    • +3
      Потому что, как говорит великодушный пожизненный диктатор, GIL, хоть и имеет слабые места, он просто не создан для того, чтобы делать многопроцессорные приложения, а скорее для работы с IO. multiprocessing для описанных в топике задач в помощь, если нужен чистый python.
      • +3
        Потоки от процессов отличаются наличием shared memory, на которую автор сделал упор. Между процессами такое придется делать руками, и желательно прямыми. Либо отказаться от общей памяти и делать общение через пайпы. Пайпы кстати автоматизированы для multiprocessing, например, очередями.
        • 0
          Ваша фантазия про отличия потоков от процессов имеет мало общего с нашей реальностью.
          Мало того, что SM изначально является механизмом IPC и к потокам отношения по сути не имеет, различия между процессами и потоками совсем не в этом.
          • 0
            Извините, я конкретно про отличия питоновских модулей Threading и Multiprocessing и их функций из коробки.
            • 0
              Ничего страшного. Но в питоновских модулях с SM тоже всё достаточно просто — multiprocessing.sharedctypes. Ну или если это называется «руками», то я уже не знаю, что значит не-руками. =)
              • 0
                Сейчас расскажу, что значит не руками. Предположим, мы отправляем в новые потоки ссылку на некий объект. Так вот в случае с потоками, в них уйдут ссылки на объект, а в случае с процессами, (дальше я точно не интересовался, как это работает, но предположу) объект запиклится и уйдет копия в каждый подпроцесс. Это я имел в виду под словом «из коробки», возможно неправильно выразился, виноват.

                Ну вот что-то типа того:
                from threading import Thread, Lock
                from multiprocessing import Process, Lock as PLock
                
                class SharedObject(object):
                    def __init__(self, property):
                        self.property = property
                
                def some_func(shared_obj, lock):
                    lock.acquire()
                    shared_obj.property += 1
                    print shared_obj.property,
                    lock.release()
                
                if __name__ == '__main__':
                    # работа с потоками:
                    lock = Lock()
                    shared_obj = SharedObject(0)
                    for thread_number in range(5):
                        thread = Thread(target=some_func, args=(shared_obj, lock))
                        thread.start()
                    # 1 2 3 4 5
                
                    # работа с процессами:
                    lock = PLock()
                    shared_obj = SharedObject(0)
                    for process_number in range(5):
                        process = Process(target=some_func, args=(shared_obj, lock))
                        process.start()
                    # 1 1 1 1 1
                
                


      • 0
        Для работы с IO GIL отпускает структуры.
        Я только не понимаю, почему нельзя использовать потокобезопасные структуры в современном мире?
        • +1
          Ох и переписывать там придется (если рассматривать CPython). Тут мне хочется надеяться что PyPy-stm будет более прогрессивным.
          • +1
            Ну так в том и суть, что python 3 это смена версии крупная, можно было и переписать, а то выходит это основной ограничитель развития языка на данный момент. Т.к. в других языках это уже решено.
            • 0
              Есть интересное видео David Beasley про анализ проблем GIL (youtube link). Там говорится, что есть надежда на улучшение GIL (не избавление, а именно исправление существующих недостатков реализации) в 3 версии Python.
  • 0
    Смотришь на этот код, и возникает ощущение, что писать на плюсах таки гораздо проще. %)
    • 0
      Автор не весь код на C++ выложил. Как минимум перегрузка * осталась за кадром.
      • –1
        Для плюсов тоже есть математические библиотеки, в которых всё нужное уже реализовано. Такого уровня велосипедостроение не нужно.
        • +9
          Я сам очень сильно люблю C++. Но вот в чем дело. Надо разделять две задачи.

          Первая задача — это промышленная реализация известного алгоритма в рамках системы с уже выработанным дизайном. И тут я бы тоже выбрал C++.

          Вторая задача — это активный поиск нужного алгоритма, его прототипирование, анализ и доработка до рабочего состояния. И здесь системы научных расчетов типа MATLAB или Scientific Python просто незаменимы. Часто за день нужно по нескольким научных статьям попробовать алгоритмы, использующие совсем разную математику. И системы расчетов предоставляют тебе готовую интерактивную среду с очень удобными абстракциями, позволяющие мговенно визуализировать результаты.
    • 0
      На плюсах писать проще, когда есть сформулированная задача, есть время продумать архитектуру, классы, вот это всё. А когда идёт научный поиск, то первоначальные абстракции могут быть совершенно неправильными, так как понимание задачи меняется в процессе решения. И вот питоне, в силу гибкости, можно довольно быстро поменять всё в программе или применить какой-нибудь dirty hack и проверить идею о которой в начале даже не задумывались (например, приватные переменные класса подкрутить). А вот когда ужё всё отлажено можно и посадить студента на C++ переписать.
  • +3
    А numexpr не пробовали? Если заменить строку в вашем numpy-варианте на R = numexpr.evaluate('1 / (1 + sqrt(Rx * Rx + Ry * Ry + Rz * Rz))'), то у меня сильно ускоряется (стаёт 0.6 секунд против 1.4 с numpy). Тут и многопоточность используется (процессор двухядерный).
    • +3
      Кстати, если то, что приведено в посте — не пример, а действительно именно такие вычисления нужны, то есть готовая реализация: scipy.spatial.distance.cdist. А именно, просто cdist(p, q) (считает просто ||p-q||) работает у меня 0.35 с, а 1/(1+cdist(p, q)) (считает искомую матрицу) — 0.6 с. Ну и, numexpr.evaluate('1/(1+d)', local_dict={'d': cdist(p, q)}) работает менее 0.5 c :)
    • 0
      Спасибо за комментарии! Numexpr действительно еще один способ параллельных вычислений на Python. Насколько я помню, numexpr использует собственную виртуальную машину, написанную на C, и собственный JIT компилятор. Как-то получилось, что я numexp редко использую. Часто было, что сложные выражение он не смог распарсить. И часто часть кода, которую хотелось распараллелить, было проще выделять в крупные блоки. Но действительно, numexpr для простых выражений очень удобный способ.

      По поводу cdist да, все верно. Но для поста я такую функцию выбрал для примера, чтобы было не слишком сложно, но вместе с тем и полезно.

      Добавлю в PS к посту эти моменты. Спасибо!
      • 0
        Если вы хотите рассмотреть всевозможные пути улучшения производительности, то попробуйте ещё Intel MKL или IPP (смотря что подходит в конкретных задачах). Использовать их из Cython очень просто, и получается весьма значительный прирост в скорости — у меня от использования только BLAS Level 1 (операции с векторами) из Intel MKL вычисления ускорились в пару раз по сравнению с наивным C/Cython кодом. Про умножение матриц даже говорить нечего.
  • +14
    Мораль: в питоне быстро работает то, что написано на C, а не на питоне.
    • –1
      Но Питон и сам написан на Си…
  • 0
    А есть результаты по pyopencl/pycuda?
    • 0
      На таких операциях (число операций линейно по объёму данных), и тем более таких размерах данных, странно ожидать прироста скорости от видеокарты. Конечно, кроме случая, когда до/после операции данные и так нужны на видеокарте.
      • 0
        Я недавно делал доклад на похожую тему — тоже с разбором полетов кода на чистом питоне, с постепенной заменой кусков на плюсы, и там же сбоку numpy для сравнения — только там задачкой было множество Мандельброта. И в числе прочих вариантов на плюсах там был вариант на C++AMP. Так вот, прирост производительности там был очень серьезный, где-то на порядок относительно многопоточной версии на PPL.
      • 0
        На днях довелось пообщаться с людьми из коммьюнити scientific Python, в т.ч. Оливье Гризелем и Брайаном Грейнджером. Так вот, в дискуссии на аналогичную тему, они сделали такое интересное замечание: производительность — это хорошо, но в рамках научных исследований важно, чтобы читатели всего этого дела могли потом легко разобраться в коде и воспроизвести те же подсчеты. Поэтому зачастую лучше использовать более медленный numpy вместо всякой экзотики, просто потому, что его знают и понимают все в данной области.
    • 0
      Для cuda, кстати, я приглядываюсь к NumbaPro (http://docs.continuum.io/numbapro/), это расширенный вариант Numba. Синтаксис там такой же простой, одной лишь аннотацией задача комплируется под CUDA:

      from numbapro import cuda
      
      @cuda.jit('void(float32[:], float32[:], float32[:])')
      def sum(a, b, result):
          i = cuda.grid(1)   # equals to threadIdx.x + blockIdx.x * blockDim.x
          result[i] = a[i] + b[i]
      

      Но тестов пока не делал.
  • НЛО прилетело и опубликовало эту надпись здесь
  • 0
    В новой версии питона появился ещё модуль statistics. Много не щупал, но то, что нужно было — работало хорошо и интуитивно по интерфейсам.
  • 0
    Интересно было бы сделать то же самое на Java (а может даже на Jython), и сравнить с результатами автора. Судя по тому, что Java хорошо умеет много потоков вокруг shared memory для потоков. Вдруг это оказалось бы эффективно, особенно если использовать готовые Java-библиотеки для расчётов (которые внутри, конечно же, сильно параллельные).
    • 0
      Для научных дел Java это ещё хуже чем C++. Да и вряд ли она будет считать быстрее чем си.
  • 0
    Раньше не слышал о Numba. Но раз уж о нем (ней?) зашла речь, то интересно было бы добавить в сравнение Nuitka. Понятно, что статью вы переписывать не будете, но буду благодарен, если в коментах поделитесь хотя бы какими-нибудь субъективными впечатлениями/ощущениями.
  • +1
    Большое спасибо за статью! Я сам физик, использую scipy/numpy в своей работе постоянно. Тоже прошёл путь MATLAB->C++/openMP->Python. Ваши примеры очень интересны, я давно хотел разобраться с cython, к сожалению, по нему очень мало документации для начинающих.

    Пробовал запустить ваши примеры, но ничего не вышло. Было бы здорово, чтобы вы написали какую-то минимальную инструкцию как это запускать. С помощью build.bat я модуль собрал, но дальше он запускаться оказывается:
    $python test_mp_python_ext.py 
    Traceback (most recent call last):
      File "test_mp_python_ext.py", line 20, in <module>
        import mp_cython_ext as mpcy
      File "-py/tests/mp_cython.pyx", line 16, in init tests.mp_cython_ext (tests/mp_cython.c:16853)
        import os, sys
    SystemError: Parent module '' not loaded, cannot perform relative import
    

    Ещё было бы хорошо, если бы вы положили в репозиторий тестовые данные, а то повторить ваши измерения невозможно. Спасибо.
    • 0
      Вы случайно не Python 3 используете? Насколько я знаю, в Python 3 были изменения в синтаксисе относительных импортов. Я пока еще использую Python 2.7.6. На Python 3 перейти пока не могу из-за библиотеки трехмерной визуализации Mayavi2, она пока сидит на Python 2.

      Тестовые данные генерируются при запуске теста на c++. Ниже краткая инструкция по запуску тестов.

      Окружение:
      Windows 8.1 x64
      Python 2.7.6
      MinGW x64-4.8.1-posix-seh-rev5

      Настройка рабочей среды:
      1. Scientific Python. Скачать и установить Anaconda 2.1.0 для Windows x64 с сайта continuum.io/downloads.
      2. MinGW compiler. Скачать MinGW x64-4.8.1-posix-seh-rev5 и установить с сайта sourceforge.net/projects/mingwbuilds/. При установке важно точно указать архитектуру x86_64 и posix threads.
      3. Добавить путь к компилятору в переменные среды. Дать команду gcc -v и убедиться, что используется именно установленный компилятор.

      Запуск на C++:
      1. Скачать ветку open-nuance с github.com/alec-kalinin/open-nuance
      2. Перейти в каталог src-cpp и в файле SConstruct в переменной CXX указать путь к установленному компилятору.
      3. Дать команду scons для сборки тестов. Если пакет scons не установлен в системе, дать команду conda install scons.
      4. Запустить benchmark-omp.exe. При этом выполнится тест на C++ и сгенерируются данные для python тестов.

      Запуск тестов:
      1. Перейти в каталог src-py и дать команду build.bat
      2. Перейти в каталог test и дать команды python test_mp_python.py, python test_mp_python_ext.py.
      • 0
        Спасибо, да у меня python3. Я попробую на втором это запустить.
  • 0
    Попробовал применить это дело в реальном проекте. Результат неутешителен. Внутри параллельного цикла openmp нельзя вызвать никакие питоновские функции. А как как сам scipy не предоставляет интерфейс cython, то никакие функции из этой библиотеки вызвать нельзя.

    Например, вам хочется посчитать значения функции Бесселя нулевого порядка в точках 1, 2, 3… и так далее.
    @cython.boundscheck(False)
    @cython.wraparound(False)
    def par(N):
        pR = np.empty(N)    
        
        cdef int i
        cdef int nN = N
        cdef double[:] R = pR
    
        with nogil, parallel():
            for i in prange(nN):
                    R[i] = scipy.special.j0(i)
                    
        return R
    


    Это не скомпилится, так как scipy.special.j0 это чисто питоновская функция. Конечно, можно вытащить Бесселя, например, из GSL и импортировать. Но зачем тогда нужна вся эта питоновская машинерия, когда уж можно на чистом C написать. Так что приведенные выше примеры работают только для простейших операций, для чего-то более серьёзного это уже не годится. Так что проблему GIL не преодолели, а замели под ковер.
    • 0
      Да, согласен. Более-менее работающая идея преодоления GIL заключается в отказе от Python объектов и в переходе на уровень С. Да, при этом действительно нельзя использовать Python методы. Но я не согласен, что это такая уж сильная проблема. Вот ведь в чем дело.

      Сами по себе NumPy и SciPy работают достаточно быстро и вполне параллельно, особенно в части линейной алгебры. Проблемы возникают тогда, когда нам приходится делать что-то нестандартное, написанное нами самостоятельно. И почти всегда в нашем самостоятельном коде можно выделить небольшую часть кода — узкое место, оптимизировав которое мы добьемся значительного прироста в производительности. Например, в моем коде это специальная процедура вычисления сингулярных интегралов.

      Очень многие качественные Python пакеты работают именно по такому принципу. Хороший пример — это пакет конечноэлементных расчетов SfePy (http://sfepy.org/) на Python. Там лишь маленькая часть всего кода оптимизирована при помощи Cython, но и этого хватает для очень даже неплохой производительности.
    • 0
      Да и вообще, а какой есть выход? На C++ экспериментировать с математикой слишком долго и непроизводительно. На MATLAB мне очень не хвататет возможностей языка и окружения, слишком все завязано на базовую линейную алгебру, шаг в сторону слишком тяжелый. SciPy идеальный вариант, но есть GIL. Но все-таки всегда можно так построить архитектуру своего приложения, что выделить небольшой низкоуровненый кусок кода и оптимизировать его при помощи Cython. И это работающий вариант.
      • 0
        Спасибо за разъяснения. Но выходов на самом деле несколько. Можно смотреть в сторону pypy, но пока numpy для него не полностью готов. Можно смотреть на новые языки типа rust и julia, там нет GIL и нативные треды. Но что меня сильно держит на python это matplotlib, уж больно приятные графики получаются без усилий.

        Сами по себе NumPy и SciPy работают достаточно быстро и вполне параллельно, особенно в части линейной алгебры.

        А как это сделать-то? Вот есть у меня СЛАУ, пишу я a=numpy.linalg.solve(M, b). Считает всё-равно в один тред, так что приходится использовать multiprocessing, что ещё те костыли, ни лямбд, ни динамических объектов.

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

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

        F(k) = \int f(x,k)dx
        


        и нам надо посчитать значения F(k) для разных k. Идеальная параллельная задача, так как f(x,k) — чистая функция. И сделать это с помощью cython не представляется возможным, потому что:
        1. Внутри nogil ни о каком scipy.integrate.quad не может быть и речи — надо писать свой интегратор или брать из GSL
        2. f(x,k) должна быть полностью реализована на cython. А это можно быть значительный кусок программы.

        Либо я не понимаю чего-то :)
        • +1
          А как это сделать-то? Вот есть у меня СЛАУ, пишу я a=numpy.linalg.solve(M, b). Считает всё-равно в один тред, так что приходится использовать multiprocessing, что ещё те костыли, ни лямбд, ни динамических объектов.

          Скорее всего, у вас Numpy не слинкован с какой-нибудь быстрой BLAS библиотекой (OpenBLAS, Intel MKL, ...) и поэтому использует свою реализацию, которая значительно медленнее (и, видимо, не параллельная). Сейчас проверил — linalg.solve работает вполне себе параллельно и грузит все ядра (у меня MKL стоит).

          2. f(x,k) должна быть полностью реализована на cython. А это можно быть значительный кусок программы.

          Если f() большая и сложная, вполне возможно, что в ней есть пара небольших кусков, которые занимают большую часть времени — их и нужно оптимизировать на Cython. Остальная же часть функции может быть также написана «на cython», но без оптимизаций и указания типов, к примеру — т.е. почти на обычном питоне.
          • 0
            Чтобы параллельно исполнять f(x,k) она должна быть объявлена как nogil, так что «почти как на обычном питоне» уже не выйдет. Тем более cython похож на python, но не равен. Особенно большие различия в классах, просто так приписать к классу cdef даже не исправляя функции не получается, особенно если есть что-то сложное типа декораторов, например.

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