3.2.4Поиск общих дуг и построение корреляционного графа

Рис. 3.4: Этап 1 - Поиск корреляций общих дуг

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

=

∑

,

(

−

)( − )

=1

∑

∑

[

( −

)]1/2[

( −

)]1/2

=1

=1

24

=							.
		∑	, = ∑
		=1			=1

∑ ∑

=1 =1

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

= ( ) Ψ Φ Θ 2,

где Ψ и Φ - количество рассматриваемых углов осевого вращения паттерна, Θ - количество рассматриваемых углов схождения сфер Эрвальда,

- количество рассматриваемых дифракционных паттернов, а - число итераций, необходимое для вычисления корреляции Пирсона, в которую неявно входит - разрешение общей дуги (число точек дуги).

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

В качестве первого подхода рассмотрим Intel Math Kernel Library. Эта библиотека была разработана специально для научных, инженерных и финансовых рассчетов. Intel Math Kernel Library (MKL) - это набор поточных и векторизованных математических действий, которые позволя-

ют ускорить всевозможные математические вычисления. Оптимизации MKL включают в себя:

∙Ускоренные операции для NumPy, SciPy, scikit-learn, и NumExpr

∙Привязка Python к низкоуровневым системным функциям MKL, что позволяет контролировать количество потоков во время исполнения

25

кода.

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

Главным образом, необходимо преобразовать все объекты np.array с которыми производятся действия в np.matrix. Библиотека MKL позволяет в автоматическом режиме применять параллельные многопоточные оптимизации матричных операций, таких как сложение, умножение, поэлементное возведение в степень и перемножение, матричное перемножение. Матричный объект np.matrix обязан иметь двумерную форму, следовательно массив, хранящий набор всех рассматриваемых профилей общих дуг отдельно взятого паттерна, и алгоритм его обработки претерпит изменения. Ранее это храшилище имело форму ( Φ × Θ × ), где Ψ и Φ - количество рассматриваемых углов осевого вращения паттерна, Θ - количество рассматриваемых углов схождения сфер Эрвальда, - разрешение общей дуги. Теперь же этот массив примет форму ( × ), где - количество рассматриваемых профилей общих дуг. Однако, поскольку в Python реализована удобная работа с подвыборками массивов, работа с этим массивом будет все так же удобной, как видно в следующем примере:

import math import numpy as np

import scipy.ndimage as ndimage

intersecton_angles = np.arange( 90., 150., 3. ) * deg2rad direction_angles = np.arange( 0., 360., 3. ) * deg2rad

#common arc resolution arc_steps = 200

#Calculate all correcation factors

def calc_weights(( i_img1, i_img2, arcs1, arcs2, )) :

# Drop the symmetry cases

if i_img2>=i_img1 : return None

26

# Enumerate for each-to-each compare

ind_grid = np.mgrid[ 0:direction_angles.size, 0:direction_angles.size ] ind = np.vstack([ ind_grid[0].reshape(-1), ind_grid[1].reshape(-1) ]).T

weights_block = np.zeros([ intersecton_angles.size, direction_angles.size * direction_angles.size ])

for i_t, theta in enumerate(intersecton_angles) :

# make matrices from arcs

m1 = np.mat( arcs1[i_t, ind[:,0]] )

m2 = np.mat( arcs2[i_t, ind[:,1]] )

#evaluating mean square difference k_as = 1. / arc_steps

weights_block[i_t] = np.array( np.sum( np.power( m1 - m2, 2. ), axis=1 ) * k_as ).reshape(-1)

#or evaluate pearson weighted factor

#weights are places by hand

w = getPearsonWeights(arc_steps) a = arcs1[i_t, ind[:,0]]

a_ = np.sum( w * a ) / np.sum( w ) b = arcs2[i_t, ind[:,0]]

b_ = np.sum( w * b ) / np.sum( w )

upper = np.sum( w * ( a - a_ ) * ( b - b_ ) )

lower = np.sqrt( np.sum( w * ( a - a_ ) ) ) * np.sqrt( np.sum( w * ( b - b_ ) ) )

weights_block[i_t] = np.array( upper/lower ).reshape(-1)

return ( i_img1, i_img2, weights_block )

Вслучае, если количество и время выполнения интерпретатора Python

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

Этот код впоследствии был протестирован на лабораторном серверном компьютере, чья конфигурация включает в себя 24 процессорных ядра. На этом компьютере был скомпилирован Python с подключением библиотеки

27

MKL. Проверив в консоли процессов автор убедился, что использование параллельных потоков имело место, однако скорость исполнения не превышала 350% по сравнению с однопоточными вычислениями. После применения упомянутого модуля multiprocessing, суммарная скорость исполнения нескольких процессов и потоков удалось привести в оптимальному уровню.

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

Рис. 3.5: Ускорение, полученное с применением MKL и multiprocessing

Возможно ли было получить еще большую оптимизацию? Короткий ответ: да, с помощью CUDA. Поскольку в распоряжении автора имелась игровая видеокарта Nvidia с поддержкой технологии вычислений CUDA, вызов на ее программирование незамедлительно принят.

Для того, чтобы реализовать в среде Python взаимодействие с видеокартой была выбрана свободная библиотека PyCUDA, которая представляет удобный интерфейс для написания шейдерных микропрограмм и их использования.

Шейдерная микропрограмма - это программа, которая вызывается вы-

28

числительными ядрами графического ускорителя и производит с ними необходимые вычисления. Графический ускоритель (GPU) - вычислительный элемент копьютера с собственной сверхвысокоскоростной оперативной памятью и большим количеством вычислительных ядер (shader processing unit). GPU предназначен для скалярных вычислений, т.е. одинаковых вычислений, производимых с большим количеством данных. Этому способствует как число процессоров, так и высокая скорость работы графической памяти, однако важно понимать, что при допущении архитектурных ошибок при написании алгоритма существенно снижается производительность вычислений. Примерами таких ошибок являются например обмен данными с видеокартой во время исполнения итераций, высокая связность процессов и их частая взаимная блокировка.

Ниже приведен пример использования библиотеки PyCUDA для вычисления Пирсоновской корреляции:

import pycuda.autoinit import pycuda.driver as drv import numpy

mod = drv.SourceModule("""

__global__ void Pearson(float *dest, float *a, float *b, float *w)

{

const int i = threadIdx.x; const int i_arc = blockIdx.x; const int amount = blockDim[0];

const int resolution = blockDim[1]; int j;

float w_sum = 0.0f; float wa_sum = 0.0f; float wb_sum = 0.0f;

for(j=0; j<amount; j++) { w_sum += w[j];

wa_sum += w[j]*a[i_arc][j]; wb_sum += w[j]*b[i_arc][j];

}

29

float a_ = wa_sum / w_sum; float b_ = wb_sum / w_sum;

float waa_sum = 0.0f; float wbb_sum = 0.0f; float waa_bb_sum = 0.0f;

for(j=0; j<amount; j++) {

waa_sum	+=	w[j]	*	(	a[i_arc][j]	-	a_	);
wbb_sum	+=	w[j]	*	(	b[i_arc][j]	-	b_	);

waa_bb_sum += w[j] * ( a[i_arc][j] - a_ ) * ( b[i_arc][j] - b_ );

}

dest[j] = waa_bb_sum * rsqrtf( waa_sum * wbb_sum )

""")

intersecton_angles = np.arange( 90., 150., 3. ) * deg2rad int_angs = intersecton_angles.size

direction_angles = np.arange( 0., 360., 3. ) * deg2rad dir_angs = direction_angles.size

# common arc resolution arc_steps = 200

pearson = mod.get_function("Pearson")

a = images[0] # look-up table for image arcs, previously mentioned as ’img_arcs_LUT’

b = images[1] # look-up table for image arcs, previously mentioned as ’img_arcs_LUT’

w = np.ones( arc_steps ) # pearson weights

dest = numpy.zeros([ int_angs * dir_angs ]) pearson(

drv.Out(dest), drv.In(a), drv.In(b), block=( int_angs * dir_angs, pearson, 1))

30

Рис. 3.6: Ускорение, полученное с применением MKL и CUDA

Добавив параллельные GPU-оптимизации в код проекта, производительность программы возросла еще в несколько десятков раз по сравнению со скоростью работы многопоточной MKL-оптимизированной программы. Ниже приведен график, отражающий производительности всех трех подходов.

Для большей наглядности график производительности представлен в линейном и логарифмическом масштабе. Как можно убедиться из приведенных данных, применение классической многопоточности вместе с MKLоптимизациями позволило достичь уменьшения скорости выполнения программы более чем на порядок (приблизительно в 20 раз). Дальнейшая работа по интеграции библиотеки PyCUDA в код проекта позволила уменьшить это число еще на два порядка, что дает суммарный прирост производительности около 2000 раз.

31

Теперь основываясь на данных о времени выполнения программы посчитаем величину производительности программы, получим следующий график:

Рис. 3.7: Производительность с применением MKL, CUDA и без них

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

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

{ 1, 1, 1} и { 2, 2, 2} достаточно будет вычислить вращение, переводящее паттерн из первой ориентации во вторую и используя значения этого относительного поворота обратиться в корреляционный массив.

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

{ , , } мы можем узнать, насколько верно данное положение паттернов друг относительно друга путем простого обращения к полученной структуре данных. В случае же, если запрашиваются углы некратные изначаль-

32

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

33