动态规划之序列联配问题

动态规划之序列联配问题

问题：Two sequence S and T,length(S)=m and length(T)=n。To identify an alignment of S and T that maximizes a scoreing function.
分析：Alignment是左右对齐的意思，经常表示产生式过程，即上面的序列S是怎么样通过下面的序列T变成的。关于序列联配问题的实际应用非常广，比如生物信息领域DNA序列的匹配问题，word里面的单词错误检测提醒等等。对于序列S和T，我们在S中添加一些空格变成S’,T中添加一些空格变成T’。我们假设有如下分原则:

d(S,T)=∑i=1nMatch(S′[i],T′[i])

这里的Match(a,b) is:
1.Match:+1, e.g.Match(‘a’, ‘a’)=1
2.Match:-1, e.g.Match(‘a’, ‘c’)=1
3.Ins/Del: -3 e.g. Match(‘C’, ‘_’)=-3
对问题总结:对于给定的S和T,问如何通过加空格使得得分最高。我们把求解过程当成一系列决策，对每个决彻部分我们决定S’[i]是怎么样通过T’[j]变化过来。我们首先考虑S的最后一个字母。它的来源有三种可能性：
one:S[m]与T[n]形成匹配，则求解原问题等于求解子问题S[1,2,3,…,m-1]与T[1,2,3,…,n-1]的对齐问题。
two:如果S[m]与空格匹配，则表示S[m]是T通过插入操作过来的，子问题班车了S[1,2,3,…,m-1]与T[1,2,3,…,n]的对齐问题。
three:如果T[n]与空格匹配，则表示S[m]是T通过删除操作过来的，子问题班车了S[1,2,3,…,m]与T[1,2,3,…,n-1]的对齐问题。
至此，我们可以得到递归表达式如下：
OPT(i,j) =

⎧⎩⎨⎪⎪Match(S[i],T[j])+OPT(i−1,j−1)Match(′−′,T[j])+OPT(i,j−1),Match(′−′,T[j])+OPT(i−1,j)

根据递归表达式我们写出源程序如下：

#include<stdio.h>#include<string.h>#include<stdlib.h>int max(int num1,int num2,int num3);int Match(char a,char b);int  Sequence_Algnment(char S[],char T[],int m,int n){    int i,j;    int (*OPT)[n+1];     OPT=(int (*)[n+1])malloc((n+1)*(m+1)*sizeof(int)); //动态申请n行m列的二维数组    for(i=0;i<m+1;i++){        OPT[i][0]=-3*i;    }    for(j=0;j<n+1;j++){        OPT[0][j]=-3*j;    }    for(i=1;i<m+1;i++){        for(j=1;j<n+1;j++){            OPT[i][j]=max( ( OPT[i-1][j-1]+Match(S[i-1],T[j-1]) ),            (OPT[i-1][j]-3),(OPT[i][j-1]-3) );        }    }    for(i=0;i<m+1;i++){        for(j=0;j<n+1;j++){            printf("%4d",OPT[i][j]);        }        printf("\n");    }    return OPT[m][n];}void main(){    int score;    char  strS[]="OCURRANCE",strT[]="OCCURRENCE";    printf("strS= %s   strT=%s\n",strS,strT);    printf("strS_length=%d\n",sizeof(strS));    printf("strT_length=%d\n",sizeof(strT));    score=Sequence_Algnment(strS,strT,9,10);    printf("score=%d \n",score);    }

虽然上述例子我们用这样的方法可以计算出来分数为4。但我们不然发现我们构造的矩阵是m by n。试想如果我们的序列都非常长，那我们开辟的m by n矩阵将是非常大。空间复杂度非常搞。对于这个问题，我们引入Hischberg算法。在计算得分矩阵的时候，我们是用初始化的列或者行去计算。如果知道了得分矩阵最左边的列，要计算下一列的分的话？我们知道任一单元的分依赖于与他临近的三个单元。所以知道第一列我们就可以算出第二列。此时，我们果断放弃第一列，计算第三列，以此类推。我们只是依赖了两个数组，就可以计算得到最后的分。实现方法如下：

int  Prefix_Space_Efficient_Sequence_Algnment (char S[],char T[],int m,int n){    int i,j;    int *score,*newscore;    score=(int *)malloc(sizeof(int)*(n+1));    newscore=(int *)malloc(sizeof(int)*(n+1));    //开辟数组    for(i=0;i<n+1;i++){        score[i]=-3*i;    }    for(i=1;i<m+1;i++){        newscore[0]=-3*i;        for(j=1;j<n+1;j++){            newscore[j]=max( (score[j-1]+Match(S[i-1],T[j-1])),             (score[j]-3) ,(newscore[j-1]-3) );        }    //  newscore[0]=0;        for(j=0;j<n+1;j++)            score[j]=newscore[j];           }    for(i=0;i<n+1;i++)        printf("%4d",score[i]);    printf("\n");    return score[n];    }

上面这个是前缀最优联配方法，后缀最优联配方法如下：

int  Suffix_Space_Efficient_Sequence_Algnment (char S[],char T[],int m,int n){    int i,j;    int *score,*newscore;    score=(int *)malloc(sizeof(int)*(m+1));    newscore=(int *)malloc(sizeof(int)*(m+1));    //开辟列数组    for(i=0;i<m+1;i++){        score[i]=-3*(m-i);    }    for(i=n-1;i>=0;i--){        newscore[m]=-3*(n-i);        for(j=m-1;j>=0;j--){            newscore[j]=max( (score[j+1]+            Match(S[i],T[j])), (score[j]-3) ,(newscore[j+1]-3) );        }        for(j=0;j<m+1;j++)            score[j]=newscore[j];           }    for(i=0;i<m+1;i++)        printf("%4d",score[i]);    printf("\n");    return score[0];    }

总结：大问题不会做，看能否整成小问题，假装小问题会做，原问题是否能做。连续决策，每一步做到最优。