12

问题：需要两个字符串之间的 LCS 长度。字符串的大小最多为 100 个字符。字母表是通常的 DNA 之一，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)（希望如此）。

然后为了寻找速度，我发现这篇文章Longest Common Subsequence 给出了Myers 的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;
}

主要有例子。例如。“tomato”和“potato”（不要评论），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 万个（如果可行的话））。

4

3 回答 3

1

我不熟悉用于计算 LCS 的比动态编程算法更高级的算法，但我想指出几点：

首先，只有在比较非常大、非常相似的字符串时，O(ND) 方法才有意义。这对你来说似乎不是这样。

其次，加速 LCD_Length 函数的渐近性能可能不是您应该关注的，因为您的字符串很短。如果您只关心寻找相似或不同的对（而不是所有对的精确 LCS），那么 Yannick 提到的 BK-tree 看起来是一种很有前途的方法。

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

于 2011-07-02T09:03:05.917 回答
1

有几种方法可以使您的计算更快：

  • 您可以使用 A* 搜索代替简单的动态编程（通过使用启发式方法，即直到 (i,j) 的部分对齐必须有 |ij| 删除或插入）。
  • 如果您将一个序列与一大堆其他序列进行比较，您可以通过计算导致该前缀的部分的动态规划表（或 A* 搜索状态）并重新使用计算的部分来节省时间. 如果您坚持使用动态编程表，您可以按字典顺序对字符串库进行排序，然后只丢弃更改的部分（例如，如果您有 'banana' 并想将其与 'panama' 和 'panamericanism' 进行比较，您可以重复使用表格中包含“panam”的部分）。
  • 如果大多数字符串非常相似，您可以通过查找公共前缀并从计算中排除公共前缀来节省时间（例如，“panama”和“panamericanism”的 LCS 是公共前缀“panam”加上“a”和“ericanism”）
  • 如果字符串非常不同，您可以使用符号计数来获得编辑次数的下限（例如，“AAAB”到“ABBB”需要至少 2 次编辑，因为一个中有 3 个 As 而只有 1 个 A另一个）。这也可以用于 A* 搜索的启发式。

编辑：对于比较相同的stings的情况，有人建议使用BK-Tree数据结构在 样本量很大时计算字符串的相似度得分的有效方式？

于 2011-07-02T08:29:35.733 回答
0

我建议获取一份 Gusfield 的关于字符串、树和序列的算法的副本，这都是关于计算生物学的字符串操作的。

于 2011-07-02T08:37:23.917 回答