dp-最长公共子序列

最长公共子序列

问题描述

最长公共子序列(LCS)是一个在一个序列集合中(通常为两个序列)用来查找所有序列中最长子序列的问题。一个数列 ,如果分别是两个或多个已知数列的子序列,且是所有符合此条件序列中最长的,则称为已知序列的最长公共子序列。
最长公共子序列问题是一个经典的计算机科学问题,也是数据比较程序,比如Diff工具,和生物信息学应用的基础。它也被广泛地应用在版本控制,比如Git用来调和文件之间的改变。

最长公共子序列不等于最长公共子串

最长公共子序列(LCS)是一个在一个序列集合中(通常为两个序列)用来查找所有序列中最长子序列的问题。这与查找最长公共子串的问题不同的地方是: 子序列不需要在原序列中占用连续的位置。 而最长公共子串(要求连续)和最长公共子序列是不同的。
另外在计算机科学中,最长递增子序列是指,在一个给定的数值序列中,找到一个子序列,使得这个子序列元素的数值依次递增,并且这个子序列的长度尽可能地大。最长递增子序列中的元素在原序列中不一定是连续的。许多与数学、算法、随机矩阵理论(英语:random matrix theory)、表示论相关的研究都会涉及最长递增子序列。解决最长递增子序列问题的算法最低要求O(n log n)的时间复杂度,这里n表示输入序列的规模。

问题分析

暴力枚举法:将两个字符串的所有子串相互比较,假设A、B两个子串的长度分别为m和n,那么两个子串分别有2^m和 2 ^n个子串,就需要比较 2 ^(n+m) ,时间复杂度特别复杂。
现在对暴力枚举做进一步改进,上面可以确定的是最长公共子序列长度一定是相同的,如果忽略空子序列的话,对于A长度为1的子序列有C(n,1)个,长度为2的子序列有C(n,2)个,……长度为n的子序列有C(n,n)个。对于B也可以做类似分析,即使只对序列A和序列B长度相同的子序列做比较,那么总的比较次数高达:

C(n,1)C(m,1)1 + C(n,2) * C(m,2) * 2+ …+C(n,p) * C(m,p)*p

其中p = min(m, n)。
在确定暴力枚举法不可取后,不妨转换下面这种思路。
我们假设两串子串(就叫A和B吧),现在假设匹配A的x位置,B的y位置(为了方便后面就叫Ax和By),那么A的(m,x)和B的(n,y)(m,n取决于两个子串的长度)已经相互匹配成功了,那现在出现的情况只有两种:匹配成功或不成功。
1、匹配成功,即Ax=By,那么最长的公共子序列长度就是在已经匹配成功的序列长度加1,LCS(x - 1, y - 1) + 1。
2、匹配不成功,即Ax!=By,这种情况下,最长的公共子序列长度不会改变并且已经求出,现在设t为最长子序列的最后一项,则t ≠ Ax和t ≠ By至少有一个成立。
2.1、如果t ≠ Ax,那么t=By,则LCS(x,y)= LCS(x – 1, y);
2.2、如果t ≠ By,那么t=Ax,则LCS(x,y) = LCS(x, y – 1)。
可是,我们事先并不知道t,由定义,我们取最大的一个,因此这种情况下,有LCS(x,y) = max(LCS(x – 1, y) , LCS(x, y – 1))。
那么可以得到递推关系:
LCS(x,y) =
(1) LCS(x - 1,y - 1) + 1 如果Ax = By
(2) max(LCS(x – 1, y) , LCS(x, y – 1)) 如果Ax != By

关于边界:一个空序列和任何序列的最长公共子序列都是空序列。

程序

// lcs:最长公共子串

#include <iostream>
#include <string.h>
#include <time.h>

using namespace std;

void solve();
void print_mat(int** len, int len1, int len2);

int main(int argc, char** argv) {
    clock_t start, end;
    start = clock();
    solve();
    end = clock();
    cout << "time cost: " << (end-start)*1000./CLOCKS_PER_SEC << " ms" << endl;

    return 0;
}

void solve() {
    char s1[] = "ACCGGTCGAGTGCGCGGAAGCCGGCCGAA";
    char s2[] = "GTCGTTCGGAATGCCGTTGCTCTGTAAA";
    int len1 = strlen(s1);
    int len2 = strlen(s2);
    char* ret = new char[min(len1, len2) + 1]{'\0'};
    int** len = new int*[len1];
    int** direct = new int*[len1];
    for (int i = 0; i < len1; ++i) {
        len[i] = new int[len2]{0};
        direct[i] = new int[len2]{0};
    }

    for (int i = 0; i < len1; ++i) {
        for (int j = 0; j < len2; ++j) {
            if (s1[i] == s2[j]) {
                len[i][j] = (i>0&&j>0) ? (len[i-1][j-1]+1) : 1;
                direct[i][j] = 3;       // 3表示向左上
            } else {
                int a = i>0 ? len[i-1][j] : 0;
                int b = j>0 ? len[i][j-1] : 0;
                if (a > b) {
                    len[i][j] = a;
                    direct[i][j] = 1;   // 1表示向上(第1维)
                } else {
                    len[i][j] = b;
                    direct[i][j] = 2;   // 2表示向左(第二维)
                }
            }
        }
    }

    cout << "len:\n";
    print_mat(len, len1, len2);
    cout << "direct:\n";
    print_mat(direct, len1, len2);

    // calc result and print
    --len1; --len2;
    while (len1 >=0 && len2 >= 0) {
        int idx = len[len1][len2];
        int dir = direct[len1][len2];
        cout << idx << ", " << dir << ", " << len1 << ", " << len2;
        switch (dir) {
            case 3:
                ret[idx-1] = s1[len1];
                cout << " " << ret[idx-1] << " ";
                --len1;
                --len2;
                break;
            case 2:
                --len2;
                break;
            case 1:
                --len1;
                break;
        }
        cout << endl;
    }

    cout << "s1:\t" << s1 << "\ns2:\t" << s2 << "\nret:\t" << ret << endl;


    for (int i = 0; i < strlen(s1); ++i) {
        delete[] len[i];
        delete[] direct[i];
    }
    delete[] len;
    delete[] direct;
    delete[] ret;
}

void print_mat(int** len, int len1, int len2) {
    for (int i = 0; i < len1; ++i) {
        for (int j = 0; j < len2; ++j) {
            cout << len[i][j] << " ";
        }
        cout << endl;
    }
}

测试:

$ g++ -o lcs lcs.cpp &7 ./lcs
len:
0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
0 0 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3
1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4
1 1 1 2 2 2 2 3 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5
1 2 2 2 3 3 3 3 4 4 4 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6
1 2 3 3 3 3 4 4 4 4 4 5 5 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7
1 2 3 4 4 4 4 5 5 5 5 5 6 6 6 7 7 7 7 7 7 7 7 8 8 8 8 8
1 2 3 4 4 4 4 5 5 6 6 6 6 6 6 7 7 7 7 7 7 7 7 8 8 9 9 9
1 2 3 4 4 4 4 5 6 6 6 6 7 7 7 7 7 7 8 8 8 8 8 8 8 9 9 9
1 2 3 4 5 5 5 5 6 6 6 7 7 7 7 7 8 8 8 8 9 9 9 9 9 9 9 9
1 2 3 4 5 5 5 6 6 6 6 7 8 8 8 8 8 8 9 9 9 9 9 10 10 10 10 10
1 2 3 4 5 5 6 6 6 6 6 7 8 9 9 9 9 9 9 10 10 10 10 10 10 10 10 10
1 2 3 4 5 5 6 7 7 7 7 7 8 9 9 10 10 10 10 10 10 10 10 11 11 11 11 11
1 2 3 4 5 5 6 7 7 7 7 7 8 9 10 10 10 10 10 11 11 11 11 11 11 11 11 11
1 2 3 4 5 5 6 7 8 8 8 8 8 9 10 11 11 11 11 11 11 11 11 12 12 12 12 12
1 2 3 4 5 5 6 7 8 8 8 8 9 9 10 11 11 11 12 12 12 12 12 12 12 12 12 12
1 2 3 4 5 5 6 7 8 9 9 9 9 9 10 11 11 11 12 12 12 12 12 12 12 13 13 13
1 2 3 4 5 5 6 7 8 9 10 10 10 10 10 11 11 11 12 12 12 12 12 12 12 13 14 14
1 2 3 4 5 5 6 7 8 9 10 10 11 11 11 11 11 11 12 12 12 12 12 13 13 13 14 14
1 2 3 4 5 5 6 7 8 9 10 10 11 12 12 12 12 12 12 13 13 13 13 13 13 13 14 14
1 2 3 4 5 5 6 7 8 9 10 10 11 12 13 13 13 13 13 13 13 14 14 14 14 14 14 14
1 2 3 4 5 5 6 7 8 9 10 10 11 12 13 14 14 14 14 14 14 14 14 15 15 15 15 15
1 2 3 4 5 5 6 7 8 9 10 10 11 12 13 14 14 14 15 15 15 15 15 15 15 15 15 15
1 2 3 4 5 5 6 7 8 9 10 10 11 12 13 14 14 14 15 16 16 16 16 16 16 16 16 16
1 2 3 4 5 5 6 7 8 9 10 10 11 12 13 14 14 14 15 16 16 17 17 17 17 17 17 17
1 2 3 4 5 5 6 7 8 9 10 10 11 12 13 14 14 14 15 16 16 17 17 18 18 18 18 18
1 2 3 4 5 5 6 7 8 9 10 10 11 12 13 14 14 14 15 16 16 17 17 18 18 19 19 19
1 2 3 4 5 5 6 7 8 9 10 10 11 12 13 14 14 14 15 16 16 17 17 18 18 19 20 20
direct:
2 2 2 2 2 2 2 2 2 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3
2 2 3 2 2 2 3 2 2 2 2 2 2 3 3 2 2 2 2 3 2 3 2 2 2 2 2 2
2 2 3 2 2 2 3 2 2 2 2 2 2 3 3 2 2 2 2 3 2 3 2 2 2 2 2 2
3 2 2 3 2 2 2 3 3 2 2 2 3 2 2 3 2 2 3 2 2 2 2 3 2 2 2 2
3 2 2 3 2 2 2 3 3 2 2 2 3 2 2 3 2 2 3 2 2 2 2 3 2 2 2 2
1 3 2 2 3 3 2 2 1 2 2 3 2 2 2 2 3 3 2 2 3 2 3 2 3 2 2 2
1 1 3 2 2 2 3 2 2 2 2 1 2 3 3 2 2 2 2 3 2 3 2 2 2 2 2 2
3 1 1 3 2 2 2 3 3 2 2 2 3 2 2 3 2 2 3 2 2 2 2 3 2 2 2 2
1 1 1 1 2 2 2 1 2 3 3 2 2 2 2 1 2 2 2 2 2 2 2 1 2 3 3 3
3 1 1 3 2 2 2 3 3 2 2 2 3 2 2 3 2 2 3 2 2 2 2 3 2 1 2 2
1 3 1 1 3 3 2 2 1 2 2 3 2 2 2 2 3 3 2 2 3 2 3 2 3 2 2 2
3 1 1 3 1 2 2 3 3 2 2 1 3 2 2 3 2 2 3 2 2 2 2 3 2 2 2 2
1 1 3 1 1 2 3 2 2 2 2 1 1 3 3 2 2 2 2 3 2 3 2 2 2 2 2 2
3 1 1 3 1 2 1 3 3 2 2 2 3 1 2 3 2 2 3 2 2 2 2 3 2 2 2 2
1 1 3 1 1 2 3 1 2 2 2 2 1 3 3 2 2 2 2 3 2 3 2 2 2 2 2 2
3 1 1 3 1 2 1 3 3 2 2 2 3 1 1 3 2 2 3 2 2 2 2 3 2 2 2 2
3 1 1 3 1 2 1 3 3 2 2 2 3 2 1 3 2 2 3 2 2 2 2 3 2 2 2 2
1 1 1 1 1 2 1 1 1 3 3 2 2 2 1 1 2 2 1 2 2 2 2 2 2 3 3 3
1 1 1 1 1 2 1 1 1 3 3 2 2 2 2 1 2 2 1 2 2 2 2 2 2 3 3 3
3 1 1 3 1 2 1 3 3 1 1 2 3 2 2 3 2 2 3 2 2 2 2 3 2 2 1 2
1 1 3 1 1 2 3 1 1 1 1 2 1 3 3 2 2 2 2 3 2 3 2 2 2 2 1 2
1 1 3 1 1 2 3 1 1 1 1 2 1 3 3 2 2 2 2 3 2 3 2 2 2 2 2 2
3 1 1 3 1 2 1 3 3 1 1 2 3 1 1 3 2 2 3 2 2 2 2 3 2 2 2 2
3 1 1 3 1 2 1 3 3 1 1 2 3 1 1 3 2 2 3 2 2 2 2 3 2 2 2 2
1 1 3 1 1 2 3 1 1 1 1 2 1 3 3 1 2 2 1 3 2 3 2 2 2 2 2 2
1 1 3 1 1 2 3 1 1 1 1 2 1 3 3 1 2 2 1 3 2 3 2 2 2 2 2 2
3 1 1 3 1 2 1 3 3 1 1 2 3 1 1 3 2 2 3 1 2 1 2 3 2 2 2 2
1 1 1 1 1 2 1 1 1 3 3 2 1 1 1 1 2 2 1 1 2 1 2 1 2 3 3 3
1 1 1 1 1 2 1 1 1 3 3 2 1 1 1 1 2 2 1 1 2 1 2 1 2 3 3 3
20, 3, 28, 27 A
19, 3, 27, 26 A
18, 2, 26, 25
18, 2, 26, 24
18, 3, 26, 23 G
17, 2, 25, 22
17, 3, 25, 21 C
16, 2, 24, 20
16, 3, 24, 19 C
15, 3, 23, 18 G
14, 2, 22, 17
14, 2, 22, 16
14, 3, 22, 15 G
13, 3, 21, 14 C
12, 3, 20, 13 C
11, 3, 19, 12 G
10, 2, 18, 11
10, 3, 18, 10 A
9, 3, 17, 9 A
8, 3, 16, 8 G
7, 3, 15, 7 G
6, 3, 14, 6 C
5, 2, 13, 5
5, 1, 13, 4
5, 1, 12, 4
5, 1, 11, 4
5, 3, 10, 4 T
4, 3, 9, 3 G
3, 1, 8, 2
3, 1, 7, 2
3, 3, 6, 2 C
2, 3, 5, 1 T
1, 3, 4, 0 G
s1:	ACCGGTCGAGTGCGCGGAAGCCGGCCGAA
s2:	GTCGTTCGGAATGCCGTTGCTCTGTAAA
ret:	GTCGTCGGAAGCCGGCCGAA
time cost: 0.335 ms
posted @ 2023-08-13 16:36  keep-minding  阅读(25)  评论(0编辑  收藏  举报