从装配线到DNA比对——神器动态规划

【前言】
对于一个问题，我们如果可以枚举所有的解，那么这个解空间我们是知道的。那么如何在解空间里面找到最优解呢？这时有一个非常好的方法，从底向上地构造整个解，而每一步都是从地层寻求最优解，这样就能保证在最终得到的一定是最优解。这就是最优子结构，有这种结构的问题，通常都可以用动态规划的办法来寻求最优解。而且它是从小规模（子问题）到大规模问题的构造，而常常这样的解法能够用一张表直观地表现出来。表中的元素是一个表达式的某个特定值，这个表达式表现的是问题和子问题的关系，也就是如何在子问题中的解中寻找最优的关系，这样的关系在例子中会非常地明了。
 
【装配流水线】
 
往往最经典的算法书里面都会讲最经典的“装配流水线”问题。因为它相对来说比较简单，一个问题，只有两个个子问题，那就是两条流水线选哪条路径。即使是最简单的例子也会有一大通的表达式搞得头有点晕，不过多看几遍应该是可以克服的。
假设一个工厂有两条装配线：装配线0 和 装配线1 ，这两个装配线中有相同数量的装配站用于装配元件。
可以用下面的图表示两个装配站：
 
 
可以用Si,j 表示第i条装配线（i 为0或1）的第j个装配站。从一个装配线转移到另一个装配线的下一站需要消耗时间Ti,j,例如从第一条线的第一站到第二条线的下一站（即第二站）所用的时间为T1,1。表明是从一线一站出发，而不用标记下一站是几，因为下一战一定是另一条线的j+1站。我们可以选择第一条线还是第二条线，来使得最终出线的时候，用时最短。
可能用语言表达很混乱，那么就用图像说明一下这些符号吧：
 点击查看原图
 
对于动态规划的解法，有一个可以遵循的套路：
  1. 理清楚是怎样获得最短用时的（最优解），这里是两条线，可以选择，要求最后用时最短
  2. 用递归来表示解的通式
  3. 计算（比较）得到最短时间（如何从两条线中选取）
  4. 回溯得到这条路径
什么也不说了，先找出最重要的递归式，然后填表：
 
F最短用时 = min(F0[n])+x1, F1[n]+x2)      x0,x1 表示最后一站到终点用时，Fi[n]表示第i条线的第n站。 比较最终谁用时最短
 
F1[j] =  min(F0[j-1]+a0,j  ,     F1[j-1] + t1,j-1 + a0,j)          表示到1线j站的最短时间是上一站是一线直接过来的，要么就是二线转移过来的，转移的时候付出的代价是t1,j-1。对于二线的时候是一样的
F1[j] =  min(F1[j-1]+a1,j  ,     F0[j-1] + t0,j-1 + a1,j)          -----------------------------1
 
而Fi[j] 比较特殊，因为一开始只能从一条路径进入，用时为ei (i 为0或1)

Fi[0] = ei + ai,0                                                                    -----------------------------2

在这里，我直接将ei也计入了第0个站中，因为进站和出站是不能换线的。

这样1式和2式，共同表达了一个递归关系，其含义也不难理解。下面用底向上方法构造一条路径，动态规划很多时候是先用某种方法求得那种是最好的，再回溯把该路径求出来。
先给出填表算法：
  
void findWay(int *a[],int *t[],int *line[],int &finalLine,int n) //n 个站点
{
	int **f = new int*[2];
	for(int k=0;k<2;k++)//总共有两条流水线，流水线0 和流水线 1
	{
		f[k] = new int[n];
	}
	f[0][0] = a[0][0];  
	f[1][0] = a[1][0];  

	for(int i=1;i<n;i++)   //从站点1开始到第n-1个站点
	{
		if(f[0][i-1]+a[0][i] <= f[1][i-1]+t[1][i-1] + a[0][i] )                //line zero	其实可以去掉共同部分a[0][i]
		{
			f[0][i] = f[0][i-1] + a[0][i];   
			line[0][i] = 0;
		}
		else
		{
			f[0][i] = f[1][i-1] + t[1][i-1] + a[0][i];
			line[0][i] = 1;
		}


		if(f[1][i-1] + a[1][i] <= f[0][i-1] + t[0][i-1] + a[1][i])			//line one
		{
			f[1][i] = f[1][i-1] + a[1][i] ;
			line[1][i] = 1;
		}
		else
		{
			f[1][i] = f[0][i-1] + t[0][i-1] + a[1][i];
			line[1][i] = 0;
		}
	}

	if(f[0][n-1] <=  f[1][n-1])
	{
		finalLine = 0; 
	}
	else
	{
		 finalLine = 1;
	}
}

然后再给出递归调用的回溯（利用填表的过程中记录的路径），得到顺序的解：

 

void printStations(int *line[],int finalLine,int n)   //递归输出经过的站点
{
	int i = finalLine;
	int j = n-1;
	if(j>0)
	printStations(line,line[i][j],j);    //递归产生顺序路径

	cout <<"装配线"<< i << "装配站"<<j<<endl;
}

 

 

假设有一个用例数据如上面的图片所示。则调用主函数

 

void main()
{
	int n = 0;
	int k = 0;			//k and g用于循环子
	int g = 0; 
	int e0 =0;			//enter
	int e1 =0;
	int x0 = 0;		//exit
	int x1 = 0;
	cout <<"总共有2条相同站点的装配线，请输入装配线的站点数 ： "<<endl;
	cin>> n;       

	cout <<"请输入进入装配线0,1的时间花费： "<<endl;
	cin>> e0 >> e1;     

	cout <<"请输入退出装配线0,1的时间花费： "<<endl;
	cin>> x0 >> x1;    
	
	int **a = new int*[2];      //站点花费时间
   int **t = new int*[2];       //换线花费时间
	int **line = new int*[2];   //记录路线轨迹

	for(k=0;k<2;k++)			//总共有两条流水线，流水线0 和流水线 1
	{
		a[k]  = new int[n];        //n个站点的流水线
		t[k] = new int[n];			//n个站点之间的交叉路线花费时间——换线时间
		line[k] = new int[n];		 //用于记录路径
	}
		

	//------------------------------------装配时间----------------------------------------
	cout <<"请输入进入站点以及在站点内装配的时间"<<endl;
	for( k=0;k<2;k++)
	{
		for(g=0;g<n;g++)
		{
			cout <<"\n输入装配线"<<k<<"的第"<<g<<"个站点所花费的时间 : ";
			cin >> a[k][g];
		}
	}
	a[0][0] += e0;			//将进入时间花费和第0号站点绑定
	a[0][n-1] += x0;		//将出线时间花费和第n-1号站点绑定 今后不用计算
	a[1][0] += e1;
	a[1][n-1] +=x1;   

	//------------------------------换线时间------------------------------------------------
   for( k=0;k<2;k++)
	{
		for(g=0;g<n-1;g++)
		{
			cout <<"\n输入装配线"<<k<<"的"<<g<<"站到另一条线的"<<g+1<<"站点所花费的时间 : ";
			cin >> t[k][g];
		}
	}
   cout << endl;

   int finalLine = 0;  //初始化
   findWay(a,t,line,finalLine,n);  //动态规划查找路径
   printStations(line,finalLine,n);  //打印路径
   system("pause");
}

 

 

 

那么结果是

 点击查看原图
【装配线的小结】
装配线问题的解空间比较小，故算法比较快就完成了求用时最少的路径。重要的是求解的过程有一个套路。那就是求递归表达式——>根据递归式填表——>回溯求得解路径。接下来就过渡到一个比较难一点的问题——DNA序列检测，就原理上来讲是万变不离其宗。
 
 
【DNA序列比对】
由于DNA序列十分庞大，故而需要非常大的内存空间来存储这个字符串，然而真正能用的软件是不可能无节制地使用内存，用来比对相对较短的序列还行。
要比对DNA的相似度最大的子串，则必须要有一种机制来评价怎样最大，因而选取一种计分的方法来评定两个字符串的相似度，越相似的得分越高。具体的评分标准是这样的，至于为什么是这样很有必要思考一下：
1、匹配一个得5分
2、不匹配一个扣1分
3、无法匹配（需要插入空格）扣掉2分   

 点击查看原图

其递归意义就是在三种情况中选最大值，max 应该在大括号前面，这里图片没有画好。

比较难理解的是无法匹配的情况，举个例子就比较好说明了
假如有一个DNA串为 GTAC
而另外一个字符串为  GAC
显然一一对应的时候（最后插入一个空格）得分是不高的，只有1分
而在第二个字符的前面插入一个空格，得到G_AC 这样得分为13分，高吧？
所以，就要考虑在哪里插入空格的复杂问题。
为了想清楚这个问题，还是用填表的方法比较直观发现问题。我们每次比对总是选取一个得分最高的方案，这个方案可能是串S1插入一个空格，或者是S2插入一个空格，或者进行比匹配了，或者进行比对不匹配。对于那些长度不想等的两个字符串，总要填入 |s1.length() - s2.length()|个空格。于是先将第0行和第0列填成 -2 * i 和 -2*j (这里i和j分别对于行和列的索引)。
 这样，每次计算要求a[i-1][j-1],a[i-1][j],a[i][j-1]这三个元素要先填好。于是我们就用从i=1，j=1 一行行填的办法可以满足这个条件。最终能够得到这样一幅图：

 点击查看原图

 

void dynamic_programming(int *a[],const string &s1,const string &s2)
{
	//initial the array 
	for(unsigned int k=0;k<=s1.length();k++)//the first colunm
	{
		a[k][0] = INDEL*k;
	}
	for(unsigned int g=0;g<=s2.length();g++)
	{
		a[0][g] = INDEL*g;
	}
	//next step: fill the table
	for(unsigned int i=1;i<=s1.length();i++)//row
	{
		for(unsigned int j=1;j<=s2.length();j++)//col
		{
			int tmp1 = 0;
			int tmp2 = 0;
			int tmp3 = 0;
			tmp1 = a[i-1][j-1] + isMatch(s1[i-1],s2[j-1]);
			tmp2 = a[i-1][j] + INDEL;
			tmp3 = a[i][j-1] + INDEL; 
			a[i][j] = max(tmp1,tmp2,tmp3);
		}
	}

	for(unsigned int r=0;r<=s1.length();r++)//row
	{
		for(unsigned int s=0;s<=s2.length();s++)//col
		{
				cout <<setw(4)<< a[r][s] ;
		}
		cout << endl;
	}
}

int max(const int &a,const int &b,const int &c)
{
	if(a > b)
	{
		if(a > c)
			return a;
		else
			return c;
	}
	else // a<b
	{
		if(b > c)
			return b;
		else
			 return c;
	}
}
int isMatch(char ch1,char ch2)
{
	if(ch1 == ch2)
		return MATCH;
	else 
		return MISMATCH;
}

 

 其中用到了一个求三个数的最大值的子函数

 

 

接下来就回溯，一一检测递归表达式中的三种情况，然后在结果中进行相应的操作.如果匹配或者不匹配，放入结果串中。如果需要插入空格，则插入空格，另外一个插入对应的S1[i]或者S2[j],因为是从对角线上找，所以一定是三种情况有且仅有一种。但是注意特殊情况，当到了i = 0 或者j=0 的时候，也就是说到了第零行或者第零列的时候，优先考虑遍历完所有两条DNA串，而不是插入空格。例如，在蓝色线上的3，虽然不匹配，但是a[1][0]那里刚好是3 + 5 = -2，造成了误解。所以我们的计分方法也是可以商榷的。这里的解决办法是INDEL（插入空格）优先。

点击查看原图

所有代码如下：

 

void trace_back(int *a[],const string &s1,const string &s2)
{
	vector <char> ans1;
	vector <char> ans2;
	ans1.push_back(' ');
	ans2.push_back(' ');
	for(int i=s1.length(),j = s2.length();i>0;)
	{
		
		 if(a[i][j] == a[i-1][j] + INDEL)
		{
			ans2.push_back('-');
			ans1.push_back(s1[i-1]);
			i--;
		}
		else if(a[i][j] == a[i][j-1] + INDEL)
		{
			ans1.push_back('-');
			ans2.push_back(s2[j-1]);
			j--;
		}
		else if(a[i][j] == a[i-1][j-1] + MATCH || a[i][j] == a[i-1][j-1]+MISMATCH) //这三个if else 语句，这条不能放到最先
		{ 
			ans1.push_back(s1[i-1]);
			ans2.push_back(s2[j-1]);
			i--;
			j--;
		}
	}

	cout << endl;
	for(unsigned int i=ans1.size()-1 ;i>0;i--)
		cout << ans1.at(i);
	cout <<endl<< endl;
	for(unsigned int j=ans2.size()-1;j>0;j--)
		cout << ans2.at(j);
	cout << endl;

	//finally print out the string
}

 

 

主函数很简单，产生两个待传入的参数即可：

 

#define INDEL -2    //insert an blank
#define MISMATCH -1  
#define MATCH  5

int main()
{
	string s2 ="AGTCCCC";
	string s1 ="ACGTATCC";
	//string s2 ="ACCGTGTCATGGGTCCACTTT";
	//string s1 ="ACAATGTTGGGTCGCTAT";
	/*string s1,s2;
	cin >> s1;
	cin >> s2;*/

	//apply a two degree matrix
	int **a = new int *[s1.length()+1]; //行			额外多一个空格
	for(unsigned int g=0;g<s1.length()+1;g++)
		a[g] = new int[s2.length()+1];	//列		额外多一个空格

	//call the function to compute the score
	dynamic_programming(a, s1, s2); 
	trace_back(a, s1, s2);

	system("pause");
	return 0;
}

 

 

如果传入的两个字符串是
//string s2 ="ACCGTGTCATGGGTCCACTTT";

//string s1 ="ACAATGTTGGGTCGCTAT";

那么可以获得这样的解：

ACAATGT--TGGGTCG-CTAT

ACCGTGTCATGGGTCCACTTT

 

 

【小结】

动态规划算法复杂度是O(N2)，从装配线到DNA序列比对，动态规划算法都表现得十分完美，总是返回给我们一个最优的解。而我们也不得不佩服算法的伟大，在现实的生活中无处不在，改变着我们的世界。动态规划算法是一项高级技术，很多书本上还有许多例子，如子串问题，矩阵连乘等问题（矩阵的运算是计算机都不能承受之重，能减少一点负担是一点）。在我看来，动态规划最重要的还是找到递归式子，找准问题和子问题的关系。找到之后就能填表了，表填好了就可以回溯得到解了。

 

by bibodeng                2012-04-21     原载于 bibodeng think beyond

posted on 2012-04-22 06:31  bibodeng  阅读(534)  评论(0)  编辑  收藏  举报

导航