文本比较算法：Needleman/Wunsch算法

时间：2021-12-09 10:53:08

本文介绍基于最长公共子序列的文本比较算法——Needleman/Wunsch算法。还是以实例说明：字符串A=kitten，字符串B=sitting那他们的最长公共子序列为ittn（注：最长公共子序列不需要连续出现，但一定是出现的顺序一致），最长公共子序列长度为4。

和LD算法类似，Needleman/Wunsch算法用的都是动态规划的思想，两者十分相似。

举例说明：A=GGATCGA，B=GAATTCAGTTA，计算LCS(A,B)。

第一步：初始化动态转移矩阵

Needleman/Wunsch算法矩阵
    G A A T T C A G T T A
  0 0 0 0 0 0 0 0 0 0 0 0
G 0                      
G 0                      
A 0                      
T 0                      
C 0                      
G 0                      
A 0                      

第二步：计算矩阵的第一行

Needleman/Wunsch算法矩阵
    G A A T T C A G T T A
  0 0 0 0 0 0 0 0 0 0 0 0
G 0 1 1 1 1 1 1 1 1 1 1 1
G 0                      
A 0                      
T 0                      
C 0                      
G 0                      
A 0                      

第三步：计算矩阵的其余各行

Needleman/Wunsch算法矩阵
    G A A T T C A G T T A
  0 0 0 0 0 0 0 0 0 0 0 0
G 0 1 1 1 1 1 1 1 1 1 1 1
G 0 1 1 1 1 1 1 1 2 2 2 2
A 0 1 2 2 2 2 2 2 2 2 2 2
T 0 1 2 2 3 3 3 3 3 3 3 3
C 0 1 2 2 3 3 4 4 4 4 4 4
G 0 1 2 2 3 3 3 4 5 5 5 5
A 0 1 2 3 3 3 3 4 5 5 5 6

则，LCS(A,B)=LCS(7,11)=6

状态转移方程是：若A(i)=B(j)，LCS(i,j)=LCS(i-1,j-1)+1；否则LCS(i,j)=max(LCS(i-1,j-1),LCS(i,j-1),LCS(i-1,j))=max(LCS(i,j-1),LCS(i-1,j))。程序实现：

/*

 *侯凯,2014-9-15

 *功能：最长子序列

 */

#include<iostream>

using namespace std;

int CalTheDistance(string A,string B)

{

    int **ptr = new int*[ A.size()+ ];

    for(int i = ; i < A.size() +  ;i++)

    {

        ptr[i] = new int[B.size() + ];

    }

    for(int i=;i<A.size()+;i++)

    {

        ptr[i][] = ;

    }

    for(int i=;i<B.size()+;i++)

    {

        ptr[][i] = ;

    }

    for(int i=;i<A.size();i++)

    {

        for(int j=;j<B.size();j++)

        {

            if(A[i]==B[j])

                ptr[i+][j+]=ptr[i][j]+;

            else

                ptr[i+][j+]=max(ptr[i+][j],ptr[i][j+]);

        }

    }

    int result = ptr[A.size()][B.size()];

    for(int i = ; i < A.size() +  ;i++)

    {

        delete [] ptr[i];

        ptr[i] = NULL;

    }

    delete[] ptr;

    ptr = NULL;

    return result;

}

int main()

{

    string str1 = "GGATCGA";

    string str2 = "GAATTCAGTTA";

    //最长子序列为6

    int distance = CalTheDistance(str1,str2);

    cout<<distance<<endl;

    system("Pause");

}

以上面为例A=GGATCGA，B=GAATTCAGTTA，LCS(A,B)=6

他们的匹配为：

A：GGA_TC_G__A

B：GAATTCAGTTA

如上面所示，蓝色表示完全匹配，黑色表示编辑操作，_表示插入字符或者是删除字符操作。如上面所示，蓝色字符有6个，表示最长公共子串长度为6。

利用上面的Needleman/Wunsch算法矩阵，通过回溯，能找到匹配字串

第一步：定位在矩阵的右下角

Needleman/Wunsch算法矩阵
    G A A T T C A G T T A
  0 0 0 0 0 0 0 0 0 0 0 0
G 0 1 1 1 1 1 1 1 1 1 1 1
G 0 1 1 1 1 1 1 1 2 2 2 2
A 0 1 2 2 2 2 2 2 2 2 2 2
T 0 1 2 2 3 3 3 3 3 3 3 3
C 0 1 2 2 3 3 4 4 4 4 4 4
G 0 1 2 2 3 3 3 4 5 5 5 5
A 0 1 2 3 3 3 3 4 5 5 5 6

第二步：回溯单元格，至矩阵的左上角

若ai=bj，则回溯到左上角单元格

Needleman/Wunsch算法矩阵
    G A A T T C A G T T A
  0 0 0 0 0 0 0 0 0 0 0 0
G 0 1 1 1 1 1 1 1 1 1 1 1
G 0 1 1 1 1 1 1 1 2 2 2 2
A 0 1 2 2 2 2 2 2 2 2 2 2
T 0 1 2 2 3 3 3 3 3 3 3 3
C 0 1 2 2 3 3 4 4 4 4 4 4
G 0 1 2 2 3 3 3 4 5 5 5 5
A 0 1 2 3 3 3 3 4 5 5 5 6

若ai≠bj，回溯到左上角、上边、左边中值最大的单元格，若有相同最大值的单元格，优先级按照左上角、上边、左边的顺序

Needleman/Wunsch算法矩阵
    G A A T T C A G T T A
  0 0 0 0 0 0 0 0 0 0 0 0
G 0 1 1 1 1 1 1 1 1 1 1 1
G 0 1 1 1 1 1 1 1 2 2 2 2
A 0 1 2 2 2 2 2 2 2 2 2 2
T 0 1 2 2 3 3 3 3 3 3 3 3
C 0 1 2 2 3 3 4 4 4 4 4 4
G 0 1 2 2 3 3 3 4 5 5 5 5
A 0 1 2 3 3 3 3 4 5 5 5 6

若当前单元格是在矩阵的第一行，则回溯至左边的单元格；若当前单元格是在矩阵的第一列，则回溯至上边的单元格

Needleman/Wunsch算法矩阵
    G A A T T C A G T T A
  0 0 0 0 0 0 0 0 0 0 0 0
G 0 1 1 1 1 1 1 1 1 1 1 1
G 0 1 1 1 1 1 1 1 2 2 2 2
A 0 1 2 2 2 2 2 2 2 2 2 2
T 0 1 2 2 3 3 3 3 3 3 3 3
C 0 1 2 2 3 3 4 4 4 4 4 4
G 0 1 2 2 3 3 3 4 5 5 5 5
A 0 1 2 3 3 3 3 4 5 5 5 6

依照上面的回溯法则，回溯到矩阵的左上角

第三步：根据回溯路径，写出匹配字串

若回溯到左上角单元格，将ai添加到匹配字串A，将bj添加到匹配字串B

若回溯到上边单元格，将ai添加到匹配字串A，将_添加到匹配字串B

若回溯到左边单元格，将_添加到匹配字串A，将bj添加到匹配字串B

搜索晚整个匹配路径，匹配字串也就完成了

可以看出，LD算法和Needleman/Wunsch算法的回溯路径是一样的。这样找到的匹配字串也是一样的。

相关文章