Параллельное программирование в Python при помощи multiprocessing и shared array

Введение.

Python замечательный язык. Связка Python + NumPy + Matplotlib, на мой взгляд, сейчас одна из лучших для научных расчётов и быстрого прототипирования алгоритмов. Но у каждого инструмента есть свои светлые и тёмные стороны. Одной из самых дискутируемых особенностей Python является GIL – Global Interpreter Lock. Я бы отнёс эту особенность к тёмной стороне инструмента. Хотя многие со мной не согласятся.

Если кратко, то GIL не позволяет в одном интерпретаторе Python эффективно использовать больше одного потока. Защитники GIL утверждают, что однопоточные программы при наличии GIL работают намного эффективнее. Но наличие GIL означает, что параллельные вычисления с использованием множества потоков и общей памяти невозможны. А это достаточно сильное ограничения в современном многоядерном мире.

Один из способов преодоления GIL при помощи потоков на C++ был недавно рассмотрен в статье: Использование Python в многопоточном приложении на C++. Я же хочу рассмотреть другой способ преодоления ограничений GIL, основанный на multiprocessing и shared array. На мой взгляд, этот способ позволяет достаточно просто и эффективно использовать процессы и разделяемую память для прозрачного параллельного программирования в стиле множества потоков и общей памяти.

Задача.

В качестве примера рассмотрим следующую задачу. В трёхмерном пространстве заданы N точек v0, v1, …, vN. Требуется для каждой пары точек вычислить функцию, зависящую от расстояния между ними. Результат будет представлять собой матрицу NxN со значениями этой функции. В качестве функции возьмем следующую: f = r^3 / 12 + r^2 / 6. Этот тест, на самом деле, не такой уж и синтетический. На вычислении таких функций от расстояния основана RBF интерполяция, которая используется во многих областях вычислительной математики.

Способ параллелизации.

В данной задаче каждая строка матрицы может вычисляться независимо. Из каждых нескольких строк матрицы сформируем независимые работы и поместим их в очередь заданий (см. рис. 1). Запустим несколько процессов. Каждый процесс будет брать из очереди следующее задание на выполнение, пока не встретит специальное задание с кодом «end». В этом случае процесс заканчивает свою работу.

Реализация на Python.

В реализации на Python у нас будут два главных метода: mpCalcDistance(nodes) и
mpCalcDistance_Worker(nodes, queue, arrD). Метод mpCalcDistance(nodes) принимает на вход список узлов, создает область общей памяти, подготавливает очередь заданий и запускает процессы. Метод mpCalcDistance_Worker(nodes, queue, arrD) это вычислительный метод, работающий в собственном потоке. Он принимает на вход список узлов, очередь заданий и область общей памяти. Рассмотрим реализацию этих методов подробнее.

Метод mpCalcDistance(nodes)

Создаем область общей памяти:

    nP = nodes.shape[0]         nQ = nodes.shape[0]      arrD = mp.RawArray(ctypes.c_double, nP * nQ) 

Создаем очередь заданий. Каждое задание не что иное, как просто диапазон строчек матрицы. Специальное значение None это признак завершения работы процесса:

    nCPU = 2     nJobs = nCPU * 36         q = nP / nJobs     r = nP % nJobs       jobs = []     firstRow = 0     for i in range(nJobs):         rowsInJob = q         if (r > 0):             rowsInJob += 1             r -= 1         jobs.append((firstRow, rowsInJob))         firstRow += rowsInJob      queue = mp.JoinableQueue()     for job in jobs:         queue.put(job)     for i in range(nCPU):         queue.put(None) 

Создаем процессы и ждем пока не закончится очередь заданий:

    workers = []     for i in range(nCPU):         worker = mp.Process(target = mpCalcDistance_Worker,                             args = (nodes, queue, arrD))         workers.append(worker)         worker.start()      queue.join() 

Из области общей памяти формируем матрицу с результатами:

    D = np.reshape(np.frombuffer(arrD), (nP, nQ))     return D 

Метод mpCalcDistance_Worker(nodes, queue, arrD)

Оборачиваем область общей памяти в матрицу:

    nP = nodes.shape[0]     nQ = nodes.shape[0]      D = np.reshape(np.frombuffer(arrD), (nP, nQ)) 

Пока не закончилась очередь заданий берем следующее задание из очереди и выполняем расчет:

    while True:         job = queue.get()         if job == None:             break          start = job[0]         stop = job[0] + job[1]                    # components of the distance vector         p = nodes[start:stop]         q = nodes.T                  Rx = p[:, 0:1] - q[0:1]         Ry = p[:, 1:2] - q[1:2]         Rz = p[:, 2:3] - q[2:3]          # calculate function of the distance         L = np.sqrt(Rx * Rx + Ry * Ry + Rz * Rz)         D[start:stop, :] = L * L * L / 12 + L * L / 6                  queue.task_done() 

Результаты

Среда выполнения: двухядерный процессор, Ubuntu 12.04, 64bit.


На первой картинке показаны времена однопоточного и двухпоточного расчета для различных N. На второй картинке показано отношение времени однопоточного расчета к двухпоточному. Видно, что начиная с N = 500 мы получаем уже существенное ускорение расчета.

Очень интересный результат получается в районе числа N = 2000. В однопоточном варианте мы получаем заметный скачок времени расчета, а при многопоточном расчете коэффициент ускорения даже превышает 2. Я объясняю это эффектом кэша. В многопоточном варианте данные для каждой задачи полностью умещаются в кэш. А в однопоточном уже нет.

Так что использование процессов и общей памяти, на мой взгляд, достаточно просто позволяет обойти ограничения GIL.

Полный код всего скрипта:

""" Python multiprocessing with shared memory example.  This example demonstrate workaround for the GIL problem. Workaround uses processes instead of threads and RawArray allocated from shared memory.  See also:     [1] http://docs.python.org/2/library/multiprocessing.html     [2] http://folk.uio.no/sturlamo/python/multiprocessing-tutorial.pdf     [3] http://www.bryceboe.com/2011/01/28/the-python-multiprocessing-queue-and-large-objects/  """  import time import ctypes import multiprocessing as mp import numpy as np import matplotlib.pyplot as plt  def generateNodes(N):     """ Generate random 3D nodes     """          return np.random.rand(N, 3)  def spCalcDistance(nodes):     """ Single process calculation of the distance function.     """          p = nodes     q = nodes.T          # components of the distance vector             Rx = p[:, 0:1] - q[0:1]     Ry = p[:, 1:2] - q[1:2]     Rz = p[:, 2:3] - q[2:3]          # calculate function of the distance     L = np.sqrt(Rx * Rx + Ry * Ry + Rz * Rz)     D = L * L * L / 12 + L * L / 6          return D      def mpCalcDistance_Worker(nodes, queue, arrD):     """ Worker process for the multiprocessing calculations     """      nP = nodes.shape[0]     nQ = nodes.shape[0]      D = np.reshape(np.frombuffer(arrD), (nP, nQ))      while True:         job = queue.get()         if job == None:             break          start = job[0]         stop = job[0] + job[1]                    # components of the distance vector         p = nodes[start:stop]         q = nodes.T                  Rx = p[:, 0:1] - q[0:1]         Ry = p[:, 1:2] - q[1:2]         Rz = p[:, 2:3] - q[2:3]          # calculate function of the distance         L = np.sqrt(Rx * Rx + Ry * Ry + Rz * Rz)         D[start:stop, :] = L * L * L / 12 + L * L / 6                  queue.task_done()     queue.task_done()  def mpCalcDistance(nodes):     """ Multiple processes calculation of the distance function.     """      # allocate shared array     nP = nodes.shape[0]         nQ = nodes.shape[0]      arrD = mp.RawArray(ctypes.c_double, nP * nQ)         # setup jobs     #nCPU = mp.cpu_count()     nCPU = 2     nJobs = nCPU * 36         q = nP / nJobs     r = nP % nJobs       jobs = []     firstRow = 0     for i in range(nJobs):         rowsInJob = q         if (r > 0):             rowsInJob += 1             r -= 1         jobs.append((firstRow, rowsInJob))         firstRow += rowsInJob      queue = mp.JoinableQueue()     for job in jobs:         queue.put(job)     for i in range(nCPU):         queue.put(None)      # run workers     workers = []     for i in range(nCPU):         worker = mp.Process(target = mpCalcDistance_Worker,                             args = (nodes, queue, arrD))         workers.append(worker)         worker.start()      queue.join()         # make array from shared memory         D = np.reshape(np.frombuffer(arrD), (nP, nQ))     return D  def compareTimes():     """ Compare execution time single processing versus multiple processing.     """     nodes = generateNodes(3000)          t0 = time.time()     spD = spCalcDistance(nodes)     t1 = time.time()     print "single process time: {:.3f} s.".format(t1 - t0)      t0 = time.time()     mpD = mpCalcDistance(nodes)     t1 = time.time()     print "multiple processes time: {:.3f} s.".format(t1 - t0)          err = np.linalg.norm(mpD - spD)     print "calculate error: {:.2e}".format(err)      def showTimePlot():         """ Generate execution time plot single processing versus multiple processing.     """          N = range(100, 4000, 4)     spTimes = []     mpTimes = []     rates = []     for i in N:         print i         nodes = generateNodes(i)                  t0 = time.time()         spD = spCalcDistance(nodes)         t1 = time.time()         sp_tt = t1 - t0         spTimes.append(sp_tt)                  t0 = time.time()         mpD = mpCalcDistance(nodes)         t1 = time.time()         mp_tt = t1 - t0         mpTimes.append(mp_tt)                  rates.append(sp_tt / mp_tt)                      plt.figure()     plt.plot(N, spTimes)     plt.plot(N, mpTimes)     plt.xlabel("N")     plt.ylabel("Execution time")          plt.figure()     plt.plot(N, rates)     plt.xlabel("N")     plt.ylabel("Rate")     plt.show()  def main():     compareTimes()     #showTimePlot()  if __name__ == '__main__':     main()