问题:需要两个字符串之间的LCS的长度。字符串的大小最多为100个字符。字母表是通常的DNA,4个字符“ACGT”。动态方法不够快。最长公共子序列长度(LCS)的快速(er)算法
我的问题是,我正在处理很多对和很多对(我可以看到数以亿计的等级)。我相信我已经将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)(希望)。
然后在寻找速度,我发现这个职位Longest Common Subsequence 这给了迈尔斯的O(ND) paper。我试过这个将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字符串是我想运行的最低字符数,最大值是100万字符(如果可行的话))。
在第二个版本中,你并没有删除'V_',如果你调用这个函数数亿次,肯定会导致问题。每次调用函数时,都应该避免分配内存,因为这会减慢你的速度。在你的情况下,我只是在栈上声明一个固定长度的数组'int V_ [2 * MAX_STRING + 2]'。 – TonyK
LCS使用的空间是O(mn)。你可以把它放到O(n),它可以让你的程序更快,因为缓存的命中率会增加。但是,如果n,m〜10^3我通常使用这个。您可以尝试我之前写过的[this](http://www.ideone.com/4rG7k)代码,并及时查看其中的差异(如果有的话)。 –
@logic_max:这是对算法的一个很好的修改,但第19行的副本将撤消避免缓存命中的任何性能优势(无论如何,我认为这对于〜40kb阵列来说并不太可怕)。您可以通过保留两个int指针来避免该副本,即在每次迭代中交换的L和L_old。 –