文本比較算法:Needleman/Wunsch算法


本文介紹基於最長公共子序列的文本比較算法——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()+ 1];
    for(int i = 0; i < A.size() + 1 ;i++)
    {
        ptr[i] = new int[B.size() + 1];
    }

    for(int i=0;i<A.size()+1;i++)
    {
        ptr[i][0] = 0;
    }
    for(int i=0;i<B.size()+1;i++)
    {
        ptr[0][i] = 0;
    }
    for(int i=0;i<A.size();i++)
    {
        for(int j=0;j<B.size();j++)
        {
            if(A[i]==B[j])
                ptr[i+1][j+1]=ptr[i][j]+1;
            else
                ptr[i+1][j+1]=max(ptr[i+1][j],ptr[i][j+1]);
        }
    }
    int result = ptr[A.size()][B.size()];
    for(int i = 0; i < A.size() + 1 ;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算法的回溯路徑是一樣的。這樣找到的匹配字串也是一樣的。


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM