动态规划（2）最优对准与最短编辑距离

 先看DNA对准

对于两个DNA序列 AACAGTTACC与TAAGGTCA

下面给出两种对准方式：

 -AACAGTTACC

TAA -GGT - - CA

AACAGTTACC

TA -AGGT- CA

当在一种对准方式中包含一个破折号（―）时，称之为插入一个缝隙。表示带有缝隙的序列存在一次缺失，或者另一序列存在一次插入。 

 前面例子中的哪种对准方式更好一点？它们都有8个匹配的碱基对上面的对准方式有两个失配碱基对，是插人了四个缝隙。而下面的对准方式中有三个失配碱基对但其代价是仅插入了两个缝隙，一般情况下如果没有事先指定失配和缝隙的罚分是不可能判定哪种对准方式更好的。例如，假定―个缝境的罚分为1，  失配的罚分为3。将一种对准方式中所有罚分之和称为该对准方式的代价。给定这些罚分分配方式：上面的对准方式的代价为10而下面的代价为1。因此上面的要好一些。但如果缝隙的罚分为2失配的罚分为1，不难看出下面的对准方式的代价较低所以也就野在为缝隙和失配指定了罚分之后就有可能确定最优对准方式了。但是如果为了做出判定而要查看所有  可能出现的对准方式．将是一件难以处理的任务。下面为序列对准问题）发一种高效的动态规划算法。

以上为《算法基础》（Richard E.Neapolitan著 第五版）的内容

对于序列x[1~n]与序列y[1~m],我们假设函数dp[i][j]表示序列x[1~i]与y[1~j]的最优对准

对于这题，很容易判断有哪几种状态通过一步操作可以得到dp[i][j]，显然是dp[i-1][j-1],dp[i-1][j],dp[i][j-1]

假设x[i]与y[j]匹配或者不匹配，那么dp[i][j]一定包含一个最优子结构dp[i-1][j-1]

假设x[i]与缝隙对准，那么dp[i][j]一定包含最优子结构dp[i-1][j](即x[1~i-1]与y[1~j]最优对准，如果x[1~i-1]与y[1~j]不是最优对准那么dp[i][j]会有更好的选择，这样与dp[i][j]的定义矛盾)

假设y[j]与缝隙对准，那么dp[i][j]一定包含最优子结构dp[i][j-1](即x[1~i]与y[1~j-1]最优对准)

因此我们可以写出dp[i][j]的递推表达式

其中x[i]与y[j]对准时的代价，当x[i]=y[j]时mismatch=0，否则则为失配罚分

gap为x[i]或y[j]与缝隙对准罚分

接下来我们要输出最后的对准序列，我们假设outputx[]与outputy[]分别是x,y对应的序列

以上分析可以知道当x[i]=y[j],并且dp[i][j]=dp[i-1][j-1]时那么对于outputx,outputy序列此时应当分别插入x[i]与y[j]

当x[i]!=y[j],并且dp[i][j]=dp[i-1][j-1]+失配罚分时那么对于outputx,outputy序列此时也应当分别插入x[i]与y[j]

当dp[i][j]=dp[i-1][j]+gap 时x[i]与缝隙对准 那么对于outputx应插入缝隙

当dp[i][j]=dp[i][j-1]+gap 时y[j]与缝隙对准 那么对于outputy应插入缝隙

最后我们要倒序输出outputx[]与outputy[]，因为我们是从表格最后一位回溯得到两个输出序列的（就像输出0-1背包中物品那样从末位回溯）

代码如下：

#include <iostream>
#include <algorithm>

#define INF 99999

using namespace std ;

int n, m;
char DNA1[101], DNA2[101];
int _mismatch, _gap ;

void alignment(      )
{
	int dp[101][101] ;//序列x1~xi与序列y1~y2最优对准

	//初始化表格
	dp[0][0] = 0 ;
	for ( int i = 1; i <= m; ++i )
		dp[0][i] = _gap*i ;//全部与缝隙对准
	for ( int i = 1; i <= n; ++i )
		dp[i][0] = _gap*i ;

	for(  int i=1; i<=n; ++i )
		for (int j = 1; j <= m; ++j) {
			int cost;
			if (DNA1[i] == DNA2[j])
				cost = 0;
			else
				cost = _mismatch ;//xi与yj失配
			int minCost = INF ;
			if ( dp[i - 1][j - 1] + cost < dp[i - 1][j] + _gap )//xi与缝隙对准
				minCost = dp[i - 1][j - 1] + cost ;
			else
				minCost = dp[i - 1][j] + _gap ;
				
			if ( dp[ i ][ j - 1 ] + _gap < minCost )
				minCost = dp[ i ][ j - 1 ] + _gap ; //yj与缝隙对准
			dp[ i ][ j ] = minCost ;
		}

	cout << dp[n][m]<<endl ;

	char outDNA1[101], outDNA2[101] ;
	int _count = 1;
	int i = n, j = m;

	while(  i>=1||j>=1  ){//注意是或不是与
			if (   DNA1[i]==DNA2[j]&& dp[i][j] == dp[i - 1][j - 1]) {//xi与yjp匹配
				outDNA1[_count] = DNA1[ i ] ;
				outDNA2[_count] = DNA2[ j ] ;
			//	cout << i << " " << j << endl;
				++_count ;
				--i, --j ;
			}
			else if (DNA1[i] != DNA2[j] && dp[i][j] == dp[i - 1][j - 1]  + _mismatch ) {
				outDNA1[_count] = DNA1[i];
				outDNA2[_count] = DNA2[j];
			//	cout << i << " " << j << endl;
				++_count;
				--i, --j;
			}
			else if ( dp[i][j] == dp[i - 1][j]+_gap ) {//xi与缝隙对准
				outDNA1[_count] = DNA1[i] ;
				outDNA2[_count] = '-' ;
				++_count;
		//		cout << i << " " << j << endl;
				--i ;
			}
			else if (dp[i][j] == dp[i][j-1] + _gap) {//yj与缝隙对准
				outDNA1[_count] = '-' ;
				outDNA2[_count] = DNA2[ j ] ;
			//	cout << i << " " << j << endl;
				++_count;
				--j ;
			}
		}

	for ( i = _count -1; i >=1; --i )//倒序输出
		cout << outDNA1[ i ] << " ";
	cout << endl ;

	for ( j = _count -1; j >=1; --j )
		cout << outDNA2[ j ] << " ";
	cout << endl ;

}


int main(     )
{
	cout << "输入两个DNA序列长度" ;
	cin >> n >> m ;

	cout << "输入两个DNA序列"<<endl;
	for (int i = 1; i <= n; ++i)
		cin >> DNA1[i] ;

	for (int i = 1; i <= m; ++i)
		cin >> DNA2[i] ;

	cout << "输入失配罚分与缝隙罚分" << endl;
	cin >> _mismatch >> _gap ;

	alignment( ) ;

	return 0 ;
}

最后的结果

下面再看最短编辑距离

编辑距离（Edit Distance），又称Levenshtein距离，是指两个字串之间，由一个转成另一个所需的最少编辑操作次数。许可的编辑操作包括将一个字符替换成另一个字符，插入一个字符，删除一个字符。一般来说，编辑距离越小，两个串的相似度越大。（以上来自百度百科）

两个问题基本相同不分析了，直接上代码，然后值得注意的是

这次输出的数据是从下标为0开始输入的，因此对于序列x[0~n-1],y[0~m-1]定义dp[i][j]为序列x[0~i-1]与y[0~j-1]的最短编辑距离。对dp[][]数组的初始化与输出需要注意一下。

#include <iostream>
#include <string>

using namespace std ;
string str1, str2;

void minedit(     )
{
	//表示x0~xi-1与y0~yi-1的最短编辑距离
	int dp[100][100];

	for (int i = 0; i <= str2.size(); ++i)
		dp[0][i] = i ;
	for (int i = 0; i <= str1.size(); ++i)
		dp[i][0] = i ;

	for(  int i=1; i<=str1.size(); ++i )
		for (int j = 1; j <= str2.size(); ++j) {
			int cost;
			if (str1[i - 1] == str2[j - 1])
				cost = 0;
			else
				cost = 1;
			int min;
			if (dp[i - 1][j - 1] + cost < dp[i - 1][j] + 1)
				min = dp[i - 1][j - 1] + cost;
			else
				min = dp[i - 1][j]+1 ;//删除xi
			if (dp[i][j - 1] + 1 < min)//添加yj
				min = dp[i][j - 1] + 1;
			dp[i][j] = min;
		}
	cout << dp[str1.size()][str2.size()] << endl ;

}

int main(   )
{
	cout << "输入两个字符串" << endl;

	while(cin >> str1>> str2 )
		minedit();


	return 0 ;
}

牛客网上有配套习题，链接：https://www.nowcoder.com/questionTerminal/9649617be3bf42288f50758df4310655

把代码略作修改即可，dp[][]大小设为dp[1100][1100]。太大太小都会发生段错误，最好还是动态吧。

最优对准参考《算法基础》（Richard E.Neapolitan著 第五版）（书中给出了分治的伪代码）