Fast (er) алгоритм для длины самой длинной общей подпоследовательности (LCS)

Проблема: нужна длина 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) документ от Myers. Я пробовал это, что связывает LCS с самым коротким редактированием script (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) также подсчитывает подстановки, но я не уверен в этом.

Любой подход, который можно реализовать (у меня нет огромных навыков программирования), приветствуется. (Если бы кто-то знал, как воспользоваться небольшим алфавитом, например).

EDIT: Я думаю, было бы полезно сказать поверх всего остального, что я использую функцию LCS между ЛЮБОЙ двумя строками в любое время. Таким образом, это не только строка say s1, по сравнению с несколькими миллионами других. Это может быть s200 с s1000, затем s0 с s10000, а затем 250 с s100000. Вероятно, не удастся также отслеживать большинство используемых строк. Я требую, чтобы длина LCS НЕ была аппроксимацией, так как я реализую алгоритм аппроксимации, и я не хочу добавлять дополнительную ошибку.

РЕДАКТОР: Просто запустил callgrind. Для ввода 10000 строк я, кажется, вызываю функцию lcs около 50 000 000 раз для разных пар строк. (10000 строк - это самый низкий, который я хочу запустить, а max - 1 миллион (если это возможно)).

Ответы

Ответ 1

Есть несколько способов ускорить вычисления:

  • Вместо простого динамического программирования вы можете использовать поиск A * (с помощью эвристики, чтобы частичное выравнивание до (i, j) обязательно должно содержать в себе /i -j | удаления или вставки).
  • Если вы сравниваете одну последовательность с целым рядом других, вы можете сэкономить время, вычислив таблицу динамического программирования (или состояние поиска A *) для части, ведущей к этому префиксу, и повторно используйте часть ваших вычислений. Если вы придерживаетесь таблицы динамического программирования, вы можете сортировать библиотеку строк по лексикографическому порядку, а затем отбрасывать только часть, которая изменяется (например, если у вас есть "банан" и вы хотите сравнить его с "панамой" и "панамериканизмом", вы можете повторно использовать часть таблицы, которая покрывает "panam" ).
  • Если большинство строк очень схожи, вы можете сэкономить время, просмотрев общий префикс и исключая общий префикс из вычисления (например, LCS из "panama" и "panamericanism" является общим префиксом "panam" plus LCS "a" и "ericanism" )
  • если строки довольно разные, вы можете использовать количество символов, чтобы получить нижнюю границу числа изменений (например, "AAAB" на "ABBB" требуется как минимум 2 редактирования, потому что есть 3 Как в одном и том же 1 A в другом). Это также можно использовать в эвристике для поиска A *.

EDIT: для случая сравнения с одним и тем же множеством один человек предложил структуру данных BK-Tree в Эффективный способ подсчета количества подобий строк при большом размере выборки?

Ответ 2

Я не знаком с алгоритмами фанатичного динамического программирования для вычислить LCS, но я хотел бы отметить несколько вещей:

Во-первых, подход O (ND) имеет смысл только в том случае, если вы сравниваете очень большие, очень аналогичные строки. Это не похоже на вас.

Во-вторых, ускорение асимптотической производительности вашей функции LCD_Length вероятно, не то, на что вы должны сосредоточиться, поскольку ваши строки довольно короткая. Если вы только заботитесь о поиске похожих или несходных пар (и не всех пары точных LCS), то BK-дерево, упомянутое Янником, выглядит многообещающим путь.

Наконец, некоторые вещи беспокоили меня о вашей реализации DP. Правильный интерпретация "board [i] [j]" в вашем коде - "самая длинная подпоследовательность строки a [1..i], b [1..j] "(я использую 1-индексирование в этих обозначениях). Поэтому, ваши петли должны включать я = lengthA и j = lengthB. Похоже, вы взломали эту ошибку в вашем коде, представив arrayLengthMacro (a), но это взлом не имеет смысла в контексте алгоритма. Когда "доска" заполнена, вы можете найти решение на плате [lengthA] [lengthB], что означает, что вы можете получить избавиться от ненужного "if (LCS < board [i] [j])" блока и возврата доска [длина Я] [lengthB]. Кроме того, границы циклов выглядят неправильно в инициализация (я уверен, что это должно быть для (i = 0; я <= maxLength; я ++) где maxLength = max (длина A, длина B)).