Проблема: Нужна длина LCS между двумя строками. Размер строк не более 100 символов. Алфавит обычный ДНК один, 4 символа "ACGT". Динамический подход недостаточно быстр.
Моя проблема в том, что я имею дело с парой лота и лота (насколько я вижу, с рангом сотен миллионов). Я считаю, что я уменьшил количество вызовов функции LCS_length до минимально возможного, поэтому единственный другой способ ускорить выполнение моей программы - использовать более эффективную функцию LCS_Length.
Я начал с реализации обычного подхода динамического программирования.
Это дает правильный ответ и, мы надеемся, будет правильно реализовано.
#define arrayLengthMacro(a) strlen(a) + 1
#define MAX_STRING 101
static int MaxLength(int lengthA, int lengthB);
/*
* Then the two strings are compared following a dynamic computing
* LCS table algorithm. Since we only require the length of the LCS
* we can get this rather easily.
*/
int LCS_Length(char *a, char *b)
{
int lengthA = arrayLengthMacro(a),lengthB = arrayLengthMacro(b),
LCS = 0, i, j, maxLength, board[MAX_STRING][MAX_STRING];
maxLength = MaxLength(lengthA, lengthB);
//printf("%d %d\n", lengthA, lengthB);
for (i = 0; i < maxLength - 1; i++)
{
board[i][0] = 0;
board[0][i] = 0;
}
for (i = 1; i < lengthA; i++)
{
for (j = 1; j < lengthB; j++)
{
/* If a match is found we allocate the number in (i-1, j-1) incremented
* by 1 to the (i, j) position
*/
if (a[i - 1] == b[j - 1])
{
board[i][j] = board[i-1][j-1] + 1;
if(LCS < board[i][j])
{
LCS++;
}
}
else
{
if (board[i-1][j] > board[i][j-1])
{
board[i][j] = board[i-1][j];
}
else
{
board[i][j] = board[i][j-1];
}
}
}
}
return LCS;
}
Это должно быть O (mn) (надеюсь).
Тогда в поисках скорости я нашел этот пост Самая длинная общая подпоследовательность
Который дал O (ND) бумагу Майерса. Я попробовал это, которое связывает LCS с кратчайшим сценарием редактирования (SES).
Соотношение, которое они дают, равно D = M + N - 2L, где D - длина SES, M и N - длины двух строк, а L - длина LCS. Но это не всегда правильно в моей реализации. Я даю свою реализацию (которую я считаю правильной):
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define arrayLengthMacro(a) strlen(a)
int LCS_Length (char *a, char *b);
int MinLength (int A, int B);
int Max (int A, int B);
int snake(int k, int max, char *a, char *b, int lengthA, int lengthB);
int main(void)
{
int L;
char a[] = "tomato", b[] = "potato"; //must give LCS = 4
L = LCS_Length(a, b);
printf("LCS: %d\n", L );
char c[] = "GTCGTTCGGAATGCCGTTGCTCTGTAAA", d[] = "ACCGGTCGAGTGCGCGGAAGCCGGCCGAA"; // must give LCS = 20
L = LCS_Length(c, d);
printf("LCS: %d\n", L );
char e[] = "human", f[] = "chimpanzee"; // must give LCS = 4
L = LCS_Length(e, f);
printf("LCS: %d\n", L );
char g[] = "howareyou", h[] = "whoareyou"; // LCS =8
L = LCS_Length(g, h);
printf("LCS: %d\n", L );
char i[] = "TTCTTTCGGTAACGCCTACTTTATGAAGAGGTTACATTGCAATCGGGTAAATTAACCAACAAGTAATGGTAGTTCCTAGTAGAGAAACCCTCCCGCTCAC",
j[] = "GCACGCGCCTGTTGCTACGCTCTGTCCATCCTTTGTGTGCCGGGTACTCAGACCGGTAACTCGAGTTGCTATCGCGGTTATCAGGATCATTTACGGACTC"; // 61
L = LCS_Length(i, j);
printf("LCS: %d\n", L );
return 0;
}
int LCS_Length(char *a, char *b)
{
int D, lengthA = arrayLengthMacro(a), lengthB = arrayLengthMacro(b),
max, *V_, *V, i, k, x, y;
max = lengthA + lengthB;
V_ = malloc(sizeof(int) * (max+1));
if(V_ == NULL)
{
fprintf(stderr, "Failed to allocate memory for LCS");
exit(1);
}
V = V_ + lengthA;
V[1] = 0;
for (D = 0; D < max; D++)
{
for (k = -D; k <= D; k = k + 2)
{
if ((k == -D && V[k-1] < V[k+1]) || (k != D && V[k-1] < V[k+1]))
{
x = V[k+1];
}
else
{
x = V[k-1] + 1;
}
y = x - k;
while ((x < lengthB) && (y < lengthA) && (a[x+1] == b[y+1]))
{
x++;
y++;
}
V[k] = x;
if ((x >= lengthB) && (y >= lengthA))
{
return (lengthA + lengthB - D)/2;
}
}
}
return (lengthA + lengthB - D)/2;
}
Есть примеры в основном.
Например. «помидор» и «картофель» (не комментируйте), длина LCS равна 4.
Реализация обнаруживает, что это 5, а SES (D в коде) получается как 2 вместо 4, что я и ожидал (удалить «t», вставить «p», удалить «m», вставить «t»). Я склонен думать, что, возможно, алгоритм O (ND) также учитывает замены, но я не уверен в этом.
Любой подход, который реализуем (у меня нет огромных навыков программирования), приветствуется.
(Если кто-то знает, как, например, воспользоваться маленьким алфавитом).
РЕДАКТИРОВАТЬ: Я думаю, что было бы полезно сказать, помимо всего прочего, что я использую функцию LCS между двумя строками в любое время. Так что это не только строка, скажем s1, по сравнению с несколькими миллионами других. Это может быть s200 с s1000, затем s0 с s10000, а затем 250 с s100000. Маловероятно, что вы сможете отслеживать наиболее часто используемые строки.
Я требую, чтобы длина LCS НЕ была аппроксимацией, поскольку я реализую алгоритм аппроксимации и не хочу добавлять дополнительную ошибку.
РЕДАКТИРОВАТЬ: Просто запустил Callgrind. Для ввода 10000 строк я, кажется, вызываю функцию lcs около 50 000 000 раз для разных пар строк. (10000 строк - это самое низкое значение, которое я хочу выполнить, и максимальное значение составляет 1 миллион (если это возможно)).