Hirschberg的LCS算法实现

解决Longest Common Subsequence(LCS)问题最常用的算法是Dyanmic programing，细节可以参考Ch15.4 of Introduction of Algorithm(2ED), MIT press, p 350。这个算法最大的问题是他的空间复杂度是O(m*n)。这样，当两个序列达到上万个节点时，内存消耗就成为了大问题。

1975年D. S. Hirschberg提出了另外一种算法，他的时间复杂度略高于Dynamic programing，但是，空间复杂度只有O(m+n)，可以很好的解决大序列的LCS问题。参见D. S. Hirschberg. A linear space algorithm for computing maximal common subsequences. Comm. A.C.M. 18(6) p341-343, 1975.

下面给出这个算法的C++和Python实现。

原算法中使用的序列下表从一开始，在此根据编程语言的特点做了优化，改成了从0开始，所以和原始算法看上去有差异。

C++(VS2005下编译通过)

 

#include  < vector >
#include  < algorithm >
using   namespace  std;

vector <   int   >  findRow0( int  m,  int  n, vector < TCHAR >  A, vector < TCHAR >  B)
... {
    vector<int> K0; 
    vector<int> K1(n+1, 0);
    //# in PDF, this lien is 1:n, It may be wrong
    for( int i = 0; i<m; i++)
    ...{
        K0 = K1;
        for(int j = 0; j < n; j++)
        ...{
            if (A[i] == B[j])
                K1[j+1] = K0[j] +1;
            else
                K1[j+1] = max ( K1[j], K0[j+1]);
        }
    }
    vector<int> LL = K1;
    
    return LL;

}

vector < TCHAR >  H_LCS0( int  m,  int  n, vector < TCHAR >   A, vector < TCHAR >  B)
... {
    
    vector<TCHAR> C;
    if (n == 0)
        C.clear(); 
    else if (m == 1)
    ...{
        vector<TCHAR>::iterator result = find( B.begin( ), B.end( ), A[0] );
        if  ( result != B.end( ) )

        //if A[0] in B:
            C = vector<TCHAR>(1,  A[0]);
        else
            C.clear(); 
    }
    else
    ...{
        int i = m / 2;
        //#step3 
        
        vector <TCHAR> A1i(A.begin(),A.begin()+i);
        vector<int> L1 = findRow0(i, n, A1i, B);
        
        
        vector <TCHAR> Anip1(A.rbegin(), A.rend()-i);
        
        vector< TCHAR > Bn1(B.rbegin(), B.rend());
        
        vector<int> L2 = findRow0(m-i, n, Anip1, Bn1);
        
        //#step4
        int M = 0;
        int k = 0;
        
        for ( int j = 0; j<=n; j++)
        ...{
            int tmp = L1[j] + L2[n-j];
            if (tmp > M)
            ...{
                M = tmp;
                k = j;
            }
        }
        
        //#step 5
        vector< TCHAR > A0i(A.begin(), A.begin()+i);
        vector< TCHAR > B0k(B.begin(), B.begin()+k);
        vector< TCHAR > C1 = H_LCS0( i, k, A0i, B0k); 
        
        vector< TCHAR > Aim(A.begin()+i, A.end());
        vector< TCHAR > Bkn(B.begin()+k, B.end());
        vector< TCHAR > C2 = H_LCS0( m-i, n-k, Aim, Bkn);


        //#step 6
        C = C1;
        C.insert(C.end(), C2.begin(), C2.end());
    }
    
    return C;
}

    

int  _tmain( int  argc, _TCHAR *  argv[])
... {
    if(argc < 3) _tprintf(_T("At least need two string "));
    else
    ...{
        int m = _tcslen(argv[1]);
        vector <TCHAR> A(argv[1], argv[1] + m);
        int n = _tcslen(argv[2]);
        vector <TCHAR> B(argv[2], argv[2] + n);
        vector <TCHAR> C = H_LCS0(m, n, A, B);
        C.push_back(0);
        _tprintf(&C[0]);

    }
    return 0;
}

Python 代码（在python2.5下测试）

def  findRow0(m, n, A, B):
     print   " findRow0 " , m , n ,  '' .join(A),  '' .join(B)
    K0  =  []
    K1  =  [0]  *  (n + 1 )
     #  in PDF, this lien is 1:n, It may be wrong
     for  i  in  range(0,m):
        K0  =  K1[:]
         for  j  in  range(0,n):
             # print i, j
             if  A[i]  ==  B[j]:
                K1[j + 1 ]  =  K0[j]  + 1
             else :
                K1[j + 1 ]  =  max ( K1[j], K0[j + 1 ])
        
    LL  =  K1
     print   ' LL = ' , LL
     return  LL

def  H_LCS0(m, n, A, B):
     print   " H_LCS0 " , m, n,  '' .join(A),  '' .join(B)
     if  n  ==  0:
        C  =  []
     elif  m  ==   1 :
         if  A[0]  in  B:
            C  =  [A[0]]
         else :
            C  =  []
     else :
        i  =  m  /   2
         # step3 
        L1  =  []
        A1i  =  A[0:i]
        L1  =  findRow0(i, n, A1i, B)
        
        
        Anip1  =  A[i:]
        Anip1.reverse()
        Bn1  =  B[:]
        Bn1.reverse()
        L2  =  findRow0(m - i, n, Anip1, Bn1)
        
         # step4
        M  =  0
        k  =  0
         for  j  in  range(0, n + 1 ):
            tmp  =  L1[j]  +  L2[n - j]
             if  tmp  >  M:
                M  =  tmp
                k  =  j
        
         # step 5
         print   ' i= ' , i,  ' k= ' , k,  ' m= ' , m,  ' n= ' , n
        C1  =  H_LCS0( i, k, A[0:i], B[0:k])
        
        C2  =  H_LCS0( m - i, n - k, A[i:], B[k:])
         # step 6
        C  =  C1 +  C2
         print   " C1= " ,   '' .join(C1),  " C2= " ,  '' .join(C2), 
     print   " C =  " ,  '' .join(C)
     return  C

A  =   " ACGTACGTACGT "
B  =   " AGTACCTACCGT "
C  =  H_LCS0(len(A), len(B), list(A), list(B))
print   " final result " ,  '' .join(C)

代码还有很多可以优化的地方。

另外，发现还有一些类似的算法，特别python的difflib采用的算法，找出的不一定是理论上的最长子序列。特别是在序列中相同元素重复出现次数比较高的时候特别明显。猜测，可能和他采用了对元素进行hash造成的。另外，他的文档中也说明：This does not yield minimal edit sequences, but does tend to yield matches that ``look right'' to people. （4.4 difflib -- Helpers for computing deltas of Python Library Reference for python 2.5）

具体算法可以参见

Pattern Matching: The Gestalt Approach
Discussion of a similar algorithm by John W. Ratcliff and D. E. Metzener. This was published in Dr. Dobb's Journal in July, 1988.