2011-07-02 85 views
10

问题:需要两个字符串之间的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万字符(如果可行的话))。

+0

在第二个版本中,你并没有删除'V_',如果你调用这个函数数亿次,肯定会导致问题。每次调用函数时,都应该避免分配内存,因为这会减慢你的速度。在你的情况下,我只是在栈上声明一个固定长度的数组'int V_ [2 * MAX_STRING + 2]'。 – TonyK

+0

LCS使用的空间是O(mn)。你可以把它放到O(n),它可以让你的程序更快,因为缓存的命中率会增加。但是,如果n,m〜10^3我通常使用这个。您可以尝试我之前写过的[this](http://www.ideone.com/4rG7k)代码,并及时查看其中的差异(如果有的话)。 –

+2

@logic_max:这是对算法的一个很好的修改,但第19行的副本将撤消避免缓存命中的任何性能优势(无论如何,我认为这对于〜40kb阵列来说并不太可怕)。您可以通过保留两个int指针来避免该副本,即在每次迭代中交换的L和L_old。 –

回答

1

有几种方法可以让你的计算速度更快:

  • 不是纯动态规划,您可以使用A *搜索(使用启发式这部分对准高达(I,J)必然有|在其中删除或插入)。
  • 如果您正在将一个序列与其他序列进行比较,则可以通过计算引导至该前缀的部分的动态编程表(或A *搜索状态)并重新使用该部分来节省时间的计算。如果你坚持使用动态编程表,你可以按字典顺序对字符串库进行排序,然后只丢弃发生变化的部分(例如,如果你有'香蕉',并想与'巴拿马'和'泛美主义'进行比较,你可以重新使用表中涵盖'panam'的部分)。
  • 如果大多数字符串非常相似,可以通过查找通用前缀并从计算中排除常用前缀来节省时间(例如,“巴拿马”和“泛美主义”的LCS是普通前缀“panam”加“a”和“ericanism”的LCS)
  • 如果字符串非常不相似,则可以使用符号计数来获得编辑次数的下限(例如,“AAAB”至“ABBB”至少需要2编辑,因为有3个在另一个中只有1个)。这也可以用于A *搜索的启发式。

编辑:为比较对的,同设置-刺的情况下,一个人建议在 Efficient way of calculating likeness scores of strings when sample size is large?

1

我不熟悉 计算LCS票友超动态编程的算法,但我想指出的几件事情:

首先,O(ND)的方法才有意义,如果你正在比较非常大,非常类似的字符串 。这似乎并不适合你。

其次,加快你的LCD_Length函数的渐近性能是 可能不是你应该关注的东西,因为你的字符串相当简短 。如果你只关心寻找类似或不相似的对(并不是所有的对都是精确的LCS),那么Yannick提到的BK树看起来像是一种很有前途的方法。

最后,有些事情让我困扰你关于DP的实现。代码中“board [i] [j]”的正确的 解释是“字符串a [1..i],b [1..j]”的最长子序列 “(我正在使用1-索引这个符号)。因此, 您的for循环应包括i = lengthA和j = lengthB。它看起来像你 通过引入arrayLengthMacro(a)在你的代码中绕过这个bug,但是在算法的上下文中这个 黑客没有意义。一旦“board”被填充,你可以在board [lengthA] [lengthB]中查找解决方案,这意味着你可以得到 摆脱不必要的“if”块并返回 板[长度A] [长度B]。此外,循环边界在初始化时看起来不对(我很确定它应该是for(i = 0; i < = maxLength; i ++) 其中maxLength = max(lengthA,lengthB))。

+0

@ Julian Panetta:谢谢你所有的观点。事实上,O(ND)不是一个很好的选择(它更像是我想的一个试验)。 第二点,我需要确切的LCS长度。 第三点,感谢您在代码中提到的内容。我有一段路要走,直到我获得编程逻辑的经验和便利。 – Yiannis

+0

噢,好吧,既然你需要所有对的LCS长度,我想你应该在比较一个字符串A和所有其他字符串B时尝试Yannick重复使用DP表的建议。构建这个的一个好方法是构建一个trie字典并为每个A运行一个DFS。然后,每次你下降trie时(从B读取一封信),你填写一行表。每次你回溯你减少你的行索引。每次你在单词树中找到一个单词节点时,就已经为(A,B)计算了一个LCS。这相当于对字符串进行排序,但使用更简洁的代码。注意:我现在使用行索引的B索引。 –

+0

我对存储东西有点犹豫,因为整个过程已经耗费空间。例如,当我有一百万个字符串时,该内存需要多少内存? – Yiannis