动态规划——序列比对问题

序列比对问题

问题要求

输入：两个序列A和B，其长度分别为m和n

输出：A和B的一种比对形式，其满足惩罚函数f(A,B)值最小。

(1) 对A和B之间的每一个空隙匹配，计惩罚分2分；

(2)对A和B之间的每一个错配，计惩罚分3分；

(3)成功配对不计惩罚分

问题分析

本题采用动态规划的思想去做

对于两个序列

A：a1,a2,a3,a4…,an

B：b1,b2,b3,b4…,bn

在两个序列中当对比到相等元素的时候自然就可以直接匹配，本题所要探究的矛盾点在于当A序列元素与B序列中的元素不等的时候到底是采取将其中的一个与空隙配对还是将两个元素进行错配，虽然空隙配对的代价会比错配的代价要小，但是如果不假思索的进行空隙配对，可能影响后面元素的配对情况

例如 ：

序列A：a c f a b c d f

序列B：a c d a b c d f

在两者进行匹配的时候如果f与d进行错配的话，难么不难看出这两个序列的惩罚函数的值是最小的，为3。但是如果采用了空隙配对的话，无论后面采用什么配对方法，惩罚函数的值较之之前的值无疑是大的。

我们知道动态规划思想实际上是分治思想的一个优化，我们先从分治的角度去理解一下本题的解题方法

假设：序列A与序列B前n-1项和前m-1项已经配对好了，要进行第an项与第bm项匹配，这时候我们考虑以下情况

an与bm相等自然直接进行匹配，惩罚函数值还是为原来的值

如果an与bm不相等是进行错配呢还是进行空隙配对呢？如果是空隙配对那么是an与空隙配对还是bm与空隙配对？（谁与空隙配对这两者有什么区别吗？）

其实是有区别的,如果an与空隙配对那么剩下的字符a1…an-1便会与b1…bn进行配对，可能这样的配对结果形成的惩罚函数会比an与bn直接错配形成的惩罚函数的值更大，反之bn与空隙配对也可能会出现最小值的情况（一下没有想好合适的例子，所以就不举例说明了）

算法设计

• 采用动态规划维护一个dp[i][j]的二维数组，该数组的含义为序列A的前i项与序列的B的前j项的惩罚函数的最小值

• dp[i][j]初始化，更具惩罚函数的定义，不难得出dp[i][0]=2*i序列A的前i项与序列的B的前0项，均为空隙匹配。同样不难得出dp[0][j]=2*j序列A的前0项与序列的B的前j项，均为空隙匹配。

• 循环遍历所有可能出现的情况（i，j的值），每次循环判断序列A的第i个与序列的B的第j个。如果两者相等，那么dp[i][j]=dp[i-1][j-1]表示二者直接匹配，惩罚函数值不累计。如果两者不等dp[i][j]=min(dp[i-1][j-1]+3,dp[i-1][j]+2,dp[i][j-1]+2)表示可能出现的情况（错配，Ai与空隙配对，Bj与空隙配对）取三者的最小值

• 最后的到dp[n][m]便为序列A与序列B的惩罚函数的最小值

程序实现

#include <iostream>
#include<cstdio>
#include<vector>
#include<stack>
using namespace std;
const int N=10005;
//string s1="abcabcab";
//string s2="abcbaca";
string s1="aabcdefgaa";
string s2="abcdefglll";
int dp[N][N];
int main()
{


    /*
       a a b c d e f g a _ a//不难看出答案是10
       _ a b c d e f g l l l

       a a b c d e f g a a _
       _ a b c d e f g l l l
    */

    int n=s1.size(),m=s2.size();//n行m列
    for(int i=0;i<=m;i++){
        dp[0][i]=2*i;
    }
    for(int i=0;i<=n;i++){
        dp[i][0]=2*i;
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=m;j++){
            if(s1[i-1]==s2[j-1])
                dp[i][j]=min(min(dp[i-1][j-1],dp[i-1][j]+2),dp[i][j-1]+2);
            else
                dp[i][j]=min(min(dp[i-1][j-1]+3,dp[i-1][j]+2),dp[i][j-1]+2);
        }
    }

    cout<<"min="<<dp[n][m]<<endl;
    return 0;
}

方案输出

从上面的程序我们不难得到结果，但是如何输出可能出现的方案呢（匹配情况），还可能有多种方案，如何输出

问题分析

• dp[i][j]的二维数组的含义为序列A的前i项与序列的B的前j项的惩罚函数的最小值，而且还有两个递推式

  dp[i][j]=dp[i-1][j-1]``（s1[i]==s2[j]）

  dp[i][j]=min(dp[i-1][j-1]+3,dp[i-1][j]+2,dp[i][j-1]+2)``（s1[i]!=s2[j]）

  通过这个递推式我们可以知道，dp[i][j]是根据前面的产生的值而生成的，选取的方式在（s1[i]==s2[j]）时有一种（s1[i]!=s2[j]）时可能会出现三种情况，我们可以根据不同的情况去判断，Ai与Bj两者此时到底是匹配的情况还是错配的情况，并且我们还可以更具不同的情况回溯到前一个位置，继续做判断直到结束

算法设计

• 声明两个动态数组v1，v2分别记录A和B可能出现的匹配序列
• 采用dfs思想对二维数组dp采用深度优先遍历，遍历出所有产生的可能情况（如果不知道dfs可以先百度学习一下）
• 采用dfs的思想遍历所有的可能，以i和j作为序列的位置，主函数开始递归dfs(n,m)
• 进入dfs，如果（s1[i]==s2[j]）表明这里是正确匹配将（s1[i]与 s2[j]）分别放入v1和v2中，递归调用dfs(i-1,j-1),回溯前一个位置，递归结束后（s1[i]与 s2[j]）分别从v1和v2中释放，因为需要遍历所有的可能。
• 如果（s1[i]！=s2[j]）依次判断三种情况，错配情况将（s1[i]与 s2[j]）分别放入v1和v2中，递归调用dfs(i-1,j-1),回溯前一个位置，递归结束后（s1[i]与 s2[j]）分别从v1和v2中释放，同样因为需要遍历所有的可能。Ai与空隙匹配情况将（s1[i]与 '_'）分别放入v1和v2中，递归调用dfs(i-1,j),回溯前一个位置(Ai与空隙匹配j的位置不用变化)，递归结束后（s1[i]与 '_'）分别从v1和v2中释放，同样因为需要遍历所有的可能。Bj与空隙匹配情况将'_' 与 s2[j]）分别放入v1和v2中，递归调用dfs(i,j-1),回溯前一个位置(Bj与空隙匹配j的位置不用变化)，递归结束后'_' 与 s2[j]）分别从v1和v2中释放，同样因为需要遍历所有的可能。
• 递归出口：i==0&j==0那么直接得出一种匹配情况，输出v1，v2；如果i==0&&j!=0表示j多了，那么剩下的Bj全部与空隙匹配，并输出v1，v2；如果i！=0&&j==0表示i多了，那么剩下的Ai全部与空隙匹配，并输出v1，v2；

源程序

#include <iostream>
#include<cstdio>
#include<vector>
#include<stack>
using namespace std;
const int N=10005;
//string s1="abcabcab";
//string s2="abcbaca";
string s1="aabcdefgaa";
string s2="abcdefglll";
int dp[N][N];
vector<char>  v1,v2;
void p(){
    for(int i=v1.size()-1;i>=0;i--)
        cout<<v1[i]<<" ";
    cout<<endl;
    for(int i=v2.size()-1;i>=0;i--)
        cout<<v2[i]<<" ";
    cout<<endl;

}
void dfs(int i,int j){
    if(i==0||j==0){
        if(i==0&&j==0)
            p();
        else if(i==0&&j!=0){
            int t=j-i;
            while(t){v1.push_back('_');v2.push_back(s2[t-1]);t--;}
            p();
            t=j-i;
            while(t--){v1.pop_back();v2.pop_back();}
        }else if(i!=0&&j==0){
            int t=i-j;
            while(t){v1.push_back(s1[t-1]);v2.push_back('_');t--;}
            p();
            t=i-j;
            while(t--){v1.pop_back();v2.pop_back();}
        }
    }
    else{
        if(s1[i-1]==s2[j-1]){
            v1.push_back(s1[i-1]);v2.push_back(s2[j-1]);
            dfs(i-1,j-1);
            v1.pop_back();v2.pop_back();
        }else{
            if(dp[i-1][j-1]+3==dp[i][j]){
                v1.push_back(s1[i-1]);v2.push_back(s2[j-1]);
                dfs(i-1,j-1);
                v1.pop_back();v2.pop_back();
            }
            if(dp[i-1][j]+2==dp[i][j]){
                v1.push_back(s1[i-1]);v2.push_back('_');
                dfs(i-1,j);
                v1.pop_back();v2.pop_back();
            }
            if(dp[i][j-1]+2==dp[i][j]){
                v1.push_back('_');v2.push_back(s2[j-1]);
                dfs(i,j-1);
                v1.pop_back();v2.pop_back();
            }
        }
    }
}
int main()
{


    /*
       a a b c d e f g a _ a
       _ a b c d e f g l l l

       a a b c d e f g a a _
       _ a b c d e f g l l l
    */

    /*


    */
    int n=s1.size(),m=s2.size();//n行m列
    for(int i=0;i<=m;i++){
        dp[0][i]=2*i;
    }
    for(int i=0;i<=n;i++){
        dp[i][0]=2*i;
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=m;j++){
            if(s1[i-1]==s2[j-1])
                dp[i][j]=min(min(dp[i-1][j-1],dp[i-1][j]+2),dp[i][j-1]+2);
            else
                dp[i][j]=min(min(dp[i-1][j-1]+3,dp[i-1][j]+2),dp[i][j-1]+2);
        }
    }
	cout<<“所有可能出现的情况”<<endl;
    dfs(n,m);
    return 0;
}

内存优化

• 我们知道，动态规划是拥有递推式的，那么前面所用的数据，对后面的作用之后那么一次，所以在只用得到值的情况之下是可以压缩存储空间的，本题中其实可以不用采取n*m大小的空间。可以直接将空间压缩为2*min(n,m)
• 可以通过两个一维数组，一个存储上一个结果的值另一个存储下一个结果的值

程序实现

#include <iostream>
#include<cstdio>
#include<vector>
#include<stack>
using namespace std;
const int N=10005;
//string s1="abcabcab";
//string s2="abcbaca";
string s1="aabcdefgaa";
string s2="abcdefglll";
int dp[N][N];
int v1[N],v2[N];
int main()
{
    int n=s1.size(),m=s2.size();//n行m列
    for(int i=0;i<=m;i++){
        v1[i]=2*i;
    }
    int i,j;
    for(int k=0;k<=m;k++){
        cout<<v1[k]<<" ";
    }
    cout<<endl;
    for(i=1;i<=n;i++){
        v2[0]=2*i;
        for(j=1;j<=m;j++){
            if(s1[i-1]==s2[j-1])
            {
                v2[j]=min(min(v1[j-1],v1[j]+2),v2[j-1]+2);
                v1[j-1]=v2[j-1];
            }else{
                v2[j]=min(min(v1[j-1]+3,v1[j]+2),v2[j-1]+2);
                v1[j-1]=v2[j-1];
            }
        }
        v1[j-1]=v2[j-1];
/*
        for(int k=0;k<=m;k++){
            cout<<v1[k]<<" ";
        }
        cout<<endl;
    }
   */
    cout<<"min="<<v1[m]<<endl;
    return 0;
}