И Снова О Gil В Python

Предисловие Область, в которой мне посчастливилось работать, называется вычислительная электрофизиология сердца .

Физиология сердечной деятельности определяется электрическими процессами, происходящими на уровне отдельных клеток миокарда.

Эти электрические процессы создают электрическое поле, которое довольно легко измерить.

Более того, оно очень хорошо описывается в рамках математических моделей электростатики.

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

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

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

Коротко о научном Python

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

Scientific Python — это набор библиотек языка Python для научных вычислений и научной визуализации.

В своей работе я использую следующие пакеты, которые покрывают примерно 90% моих потребностей:

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

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

Я не буду касаться вопросов использования CUDA и систем с раздельной памятью (чаще всего используется стандарт MPI).

Проблема в GIL. GIL (Global Interpreter Lock) — это блокировка (мьютекс), которая не позволяет нескольким потокам выполнять один и тот же байт-код. К сожалению, эта блокировка необходима, поскольку система управления памятью CPython не является потокобезопасной.

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

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

Тестовое задание

Необходимо построить матрицу р измерение Н х Н , каждый элемент р _{я, Джей} который рассчитывается по формуле:



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

Эта матрица довольно часто используется в реальных расчетах, например, при RBF-интерполяции или решении дифуров в PE методом интегральных уравнений.

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

Результаты основаны на среднем времени 10 пробежек.

Полную реализацию тестовых задач можно найти на GitHub [2].

Правильная точка в комментариях от "@chersaya".

Эта тестовая задача используется здесь в качестве примера.

Если вам действительно нужно вычислить попарные расстояния, лучше использовать функцию scipy.spatial.distance.cdist.

Параллельная реализация на 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));
         }
     }
 }
 

   

Перегруженный оператор «*» в этом классе имеет значение скалярного произведения.

Для представления набора векторов я использовал std::vector. Технология OpenMP использовалась для параллельных вычислений.

Для распараллеливания алгоритма достаточно использовать директиву #pragma omp Parallel for. Полученные результаты:

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

1.Наивная реализация на чистом 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

Здесь п , д — входные данные в формате NumPy с массивами измерений (Н, 3) И (3, Н) .

А дальше идет честный цикл на Python, вычисляющий элементы матрицы R. Полученные результаты:

В общем, это совсем не вариант численных расчетов.

2 Однопроцессорный 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. Полученные результаты:

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

Но пока это все результаты одного процессора.

Перейдем к многопроцессорности.

3 Мультипроцессор NumPy

Все бы ничего, но есть проблема.

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

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

Я точно не знаю, что стоит за RawArray. Мне кажется, что это файлы, отображаемые в памяти.

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

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

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

Полученные результаты:

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

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

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

Но это не так страшно, как кажется.

4 Китон

Модули 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. На мой взгляд, это довольно просто и красиво! Полученные результаты:

Отставание примерно в 1,1 раза как в однопроцессорном, так и в многопроцессорном варианте практически незаметно в реальных задачах.

5 Нумба

Идея здесь примерно такая же, как и в 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

Мультипроцессор Нумба:

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 была попытка реализовать параллельный цикл с помощью одной инструкции, но из-за проблем со стабильностью эта возможность была удалена в последующих выпусках.

Я уверен, что эта особенность будет исправлена со временем.

Результаты выполнения:

выводы

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



Рис.

2. Результаты многопроцессорных расчетов Мне кажется, что проблемы GIL в Python для численных расчетов практически преодолены.

На данный момент я бы рекомендовал Cython в качестве технологии параллельных вычислений.

Но я бы очень внимательно присмотрелся к Нумбе.

Ссылки

Это использует библиотеку 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)")

Теги: #научный Python #gil #параллельное программирование #python #Параллельное программирование

• Грызуны

  19 Oct, 24

• Cognitive Technologies Запустила Исследование Моральности Беспилотного Автомобиля В Случае Аварии

  19 Oct, 24

• Выпускаем Автомобиль В Интернет! Модернизация Радиостанции Prology Mdn-2740T И Аналоги

  19 Oct, 24

• Как Пережить Ядерный Апокалипсис

  19 Oct, 24

• У Вас Дома Есть Кот?

  19 Oct, 24

• Создание Форума На Drupal

  19 Oct, 24

• «Вредные» Клиенты

  19 Oct, 24

• Система Мотивации Разработчиков

  19 Oct, 24

• Команда Lingualeo Едет В Черногорию На 2 Месяца. На Борту 5 Свободных Мест! Миссия Команды — Выпускать Самые Важные Релизы

  19 Oct, 24

• 3D-Соцсети Frenzoo Понравились Основателям Skype

  19 Oct, 24