Как ускорить расчет расстояния Левенштейн
Я пытаюсь запустить симуляцию для проверки среднего расстояния Левенштейна между случайными
двоичные строки.
Моя программа находится в python, но я использую этот расширение C. Функция, которая имеет значение и занимает большую часть времени, вычисляет расстояние Левенштейна между двумя строками и является ли это.
lev_edit_distance(size_t len1, const lev_byte *string1,
size_t len2, const lev_byte *string2,
int xcost)
{
size_t i;
size_t *row; /* we only need to keep one row of costs */
size_t *end;
size_t half;
/* strip common prefix */
while (len1 > 0 && len2 > 0 && *string1 == *string2) {
len1--;
len2--;
string1++;
string2++;
}
/* strip common suffix */
while (len1 > 0 && len2 > 0 && string1[len1-1] == string2[len2-1]) {
len1--;
len2--;
}
/* catch trivial cases */
if (len1 == 0)
return len2;
if (len2 == 0)
return len1;
/* make the inner cycle (i.e. string2) the longer one */
if (len1 > len2) {
size_t nx = len1;
const lev_byte *sx = string1;
len1 = len2;
len2 = nx;
string1 = string2;
string2 = sx;
}
/* check len1 == 1 separately */
if (len1 == 1) {
if (xcost)
return len2 + 1 - 2*(memchr(string2, *string1, len2) != NULL);
else
return len2 - (memchr(string2, *string1, len2) != NULL);
}
len1++;
len2++;
half = len1 >> 1;
/* initalize first row */
row = (size_t*)malloc(len2*sizeof(size_t));
if (!row)
return (size_t)(-1);
end = row + len2 - 1;
for (i = 0; i < len2 - (xcost ? 0 : half); i++)
row[i] = i;
/* go through the matrix and compute the costs. yes, this is an extremely
* obfuscated version, but also extremely memory-conservative and relatively
* fast. */
if (xcost) {
for (i = 1; i < len1; i++) {
size_t *p = row + 1;
const lev_byte char1 = string1[i - 1];
const lev_byte *char2p = string2;
size_t D = i;
size_t x = i;
while (p <= end) {
if (char1 == *(char2p++))
x = --D;
else
x++;
D = *p;
D++;
if (x > D)
x = D;
*(p++) = x;
}
}
}
else {
/* in this case we don't have to scan two corner triangles (of size len1/2)
* in the matrix because no best path can go throught them. note this
* breaks when len1 == len2 == 2 so the memchr() special case above is
* necessary */
row[0] = len1 - half - 1;
for (i = 1; i < len1; i++) {
size_t *p;
const lev_byte char1 = string1[i - 1];
const lev_byte *char2p;
size_t D, x;
/* skip the upper triangle */
if (i >= len1 - half) {
size_t offset = i - (len1 - half);
size_t c3;
char2p = string2 + offset;
p = row + offset;
c3 = *(p++) + (char1 != *(char2p++));
x = *p;
x++;
D = x;
if (x > c3)
x = c3;
*(p++) = x;
}
else {
p = row + 1;
char2p = string2;
D = x = i;
}
/* skip the lower triangle */
if (i <= half + 1)
end = row + len2 + i - half - 2;
/* main */
while (p <= end) {
size_t c3 = --D + (char1 != *(char2p++));
x++;
if (x > c3)
x = c3;
D = *p;
D++;
if (x > D)
x = D;
*(p++) = x;
}
/* lower triangle sentinel */
if (i <= half) {
size_t c3 = --D + (char1 != *char2p);
x++;
if (x > c3)
x = c3;
*p = x;
}
}
}
i = *end;
free(row);
return i;
}
Можно ли это ускорить?
Я буду запускать код в 32-разрядном ubuntu на восьмипроцессорном процессоре AMD FX (tm) -8350.
Вот код python, который вызывает его.
from Levenshtein import distance
import random
for i in xrange(16):
sum = 0
for j in xrange(1000):
str1 = bin(random.getrandbits(2**i))[2:].zfill(2**i)
str2 = bin(random.getrandbits(2**i))[2:].zfill(2**i)
sum += distance(str1,str2)
print i,sum/(1000*2**i)
Ответы
Ответ 1
Возможно, вы могли бы запустить эту параллель. Сгенерируйте один гигантский список рандомов в начале, затем в вашем цикле, создайте потоки (8 потоков) за один раз для каждого процесса одним фрагментом списка и добавьте его окончательный результат в переменную суммы. Или сгенерируйте список из 8 одновременно и сделайте 8 за раз.
Проблема с предложением openmp: "Этот алгоритм плохо распараллеливается из-за большого количества зависимостей данных" - Wikipedia
from threading import Thread
sum = 0
def calc_distance(offset) :
sum += distance(randoms[offset][0], randoms[offset][1]) #use whatever addressing scheme is best
threads = []
for i in xrange(8) :
t = new Thread(target=calc_distance, args=(i))
t.start()
threads.append(t)
позже....
for t in threads :
t.join()
Я думаю, что этот метод будет хорошо переноситься на opencl позже, если бы ядро levenshtein было доступно (или кодируемое).
Это просто быстрый пост из памяти, поэтому, возможно, есть некоторые изгибы.
Ответ 2
Что вы можете сделать, это начать с изучения некоторых концепций и директив OpenMP с этого сайта: Начальный Primer для OpenMP
Вам нужен компилятор, совместимый с OpenMP. Вот список компиляторов, который работает. При компиляции кода вы захотите использовать параметр -fopenmp
.
Я только добавил директиву компилятора #pragma omp parallel for
к вашему коду, чтобы сообщить компилятору, что следующие блоки кода могут выполняться параллельно. Вы можете увидеть прирост производительности в результате изменения циклов while и циклов или применения шаблона OpenMP во всей этой функции. Вы можете настроить производительность, отрегулировав количество потоков, которые используются для выполнения циклов for, используя функцию omp_set_num_threads()
перед этими блоками. Хорошее количество для вас - это 8, так как вы будете работать на 8-ядерном процессоре.
lev_edit_distance(size_t len1, const lev_byte *string1,
size_t len2, const lev_byte *string2,
int xcost)
{
size_t i;
size_t *row; /* we only need to keep one row of costs */
size_t *end;
size_t half;
// Set the number of threads the OpenMP framework will use to parallelize the for loops
omp_set_num_threads(8);
/* strip common prefix */
while (len1 > 0 && len2 > 0 && *string1 == *string2) {
len1--;
len2--;
string1++;
string2++;
}
/* strip common suffix */
while (len1 > 0 && len2 > 0 && string1[len1-1] == string2[len2-1]) {
len1--;
len2--;
}
/* catch trivial cases */
if (len1 == 0)
return len2;
if (len2 == 0)
return len1;
/* make the inner cycle (i.e. string2) the longer one */
if (len1 > len2) {
size_t nx = len1;
const lev_byte *sx = string1;
len1 = len2;
len2 = nx;
string1 = string2;
string2 = sx;
}
/* check len1 == 1 separately */
if (len1 == 1) {
if (xcost)
return len2 + 1 - 2*(memchr(string2, *string1, len2) != NULL);
else
return len2 - (memchr(string2, *string1, len2) != NULL);
}
len1++;
len2++;
half = len1 >> 1;
/* initalize first row */
row = (size_t*)malloc(len2*sizeof(size_t));
if (!row)
return (size_t)(-1);
end = row + len2 - 1;
#pragma omp parallel for
for (i = 0; i < len2 - (xcost ? 0 : half); i++)
row[i] = i;
/* go through the matrix and compute the costs. yes, this is an extremely
* obfuscated version, but also extremely memory-conservative and relatively
* fast. */
if (xcost) {
#pragma omp parallel for
for (i = 1; i < len1; i++) {
size_t *p = row + 1;
const lev_byte char1 = string1[i - 1];
const lev_byte *char2p = string2;
size_t D = i;
size_t x = i;
while (p <= end) {
if (char1 == *(char2p++))
x = --D;
else
x++;
D = *p;
D++;
if (x > D)
x = D;
*(p++) = x;
}
}
}
else {
/* in this case we don't have to scan two corner triangles (of size len1/2)
* in the matrix because no best path can go throught them. note this
* breaks when len1 == len2 == 2 so the memchr() special case above is
* necessary */
row[0] = len1 - half - 1;
#pragma omp parallel for
for (i = 1; i < len1; i++) {
size_t *p;
const lev_byte char1 = string1[i - 1];
const lev_byte *char2p;
size_t D, x;
/* skip the upper triangle */
if (i >= len1 - half) {
size_t offset = i - (len1 - half);
size_t c3;
char2p = string2 + offset;
p = row + offset;
c3 = *(p++) + (char1 != *(char2p++));
x = *p;
x++;
D = x;
if (x > c3)
x = c3;
*(p++) = x;
}
else {
p = row + 1;
char2p = string2;
D = x = i;
}
/* skip the lower triangle */
if (i <= half + 1)
end = row + len2 + i - half - 2;
/* main */
while (p <= end) {
size_t c3 = --D + (char1 != *(char2p++));
x++;
if (x > c3)
x = c3;
D = *p;
D++;
if (x > D)
x = D;
*(p++) = x;
}
/* lower triangle sentinel */
if (i <= half) {
size_t c3 = --D + (char1 != *char2p);
x++;
if (x > c3)
x = c3;
*p = x;
}
}
}
i = *end;
free(row);
return i;
}
Вы также можете выполнять операции reduction переменных, которые также используются в ваших циклах for, чтобы обеспечить простые параллельные вычисления, например суммирование, умножение и т.д.
int main()
{
int i = 0,
j = 0,
sum = 0;
char str1[30]; // Change size to fit your specifications
char str2[30];
#pragma omp parallel for
for(i=0;i<16;i++)
{
sum = 0;
// Could do a reduction on sum across all threads
for(j=0;j<1000;j++)
{
// Calls will have to be changed
// I don't know much Python so I'll leave that to the experts
str1 = bin(random.getrandbits(2**i))[2:].zfill(2**i)
str2 = bin(random.getrandbits(2**i))[2:].zfill(2**i)
sum += distance(str1,str2)
}
printf("%d %d",i,(sum/(1000*2*i)));
}
}
Ответ 3
Что я буду делать:
1) Очень небольшая оптимизация: выделяйте раз и навсегда row
, чтобы избежать накладных расходов на управление памятью. Или вы можете попробовать realloc()
, или вы можете отслеживать размер row
в статической переменной (и также иметь статику row
). Тем не менее, это экономит очень мало, даже если это будет немного дешевле.
2) Вы пытаетесь вычислить среднее значение. Сделайте средний расчет и на C. Это должно что-то сэкономить в звонках. Опять же, небольшие изменения, но это дешево.
3) Поскольку вас не интересуют фактические вычисления, но только в результатах, то, скажем, у вас есть три компьютера, и каждый из них является четырехъядерным. Затем запустите на каждом из них четыре экземпляра программы, причем цикл будет в 12 раз короче. В двенадцатый раз вы получите двенадцать результатов: в среднем это, а Боб - ваш дядя.
Вариант № 3 не требует никаких изменений, кроме цикла, и вы можете сделать его параметром командной строки, чтобы вы могли развернуть программу на переменном количестве компьютеров. Фактически, вы можете вывести как результат, так и его "вес", чтобы минимизировать вероятность ошибок при суммировании результатов.
for j in xrange(N):
str1 = bin(random.getrandbits(2**i))[2:].zfill(2**i)
str2 = bin(random.getrandbits(2**i))[2:].zfill(2**i)
sum += distance(str1,str2)
print N,i,sum/(N*2**i)
Но если вас интересует общая статистика Levenshtein, я не уверен, что выполнение вычислений с использованием только 0 и 1 символов подходит для вашей цели. Из строки 01010101 вы получаете 10101010 либо путем переключения восьми символов, либо путем сброса первого и добавления нуля в конце с двумя разными расходами. Если у вас есть все буквы алфавита, вторая возможность становится гораздо менее вероятной, и это должно что-то изменить в средневзвешенном сценарии. Или я что-то упускаю?
Ответ 4
Кто-то еще провел много исследований год или два назад и также проверил время выполнения.
Он придумал этот и в основном использовал дерево решений, чтобы ускорить процесс.