从装配线到DNA比对——神器动态规划
- 理清楚是怎样获得最短用时的(最优解),这里是两条线,可以选择,要求最后用时最短
- 用递归来表示解的通式
- 计算(比较)得到最短时间(如何从两条线中选取)
- 回溯得到这条路径
Fi[0] = ei + ai,0 -----------------------------2
在这里,我直接将ei也计入了第0个站中,因为进站和出站是不能换线的。
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"); }
那么结果是
其递归意义就是在三种情况中选最大值,max 应该在大括号前面,这里图片没有画好。
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序列比对,动态规划算法都表现得十分完美,总是返回给我们一个最优的解。而我们也不得不佩服算法的伟大,在现实的生活中无处不在,改变着我们的世界。动态规划算法是一项高级技术,很多书本上还有许多例子,如子串问题,矩阵连乘等问题(矩阵的运算是计算机都不能承受之重,能减少一点负担是一点)。在我看来,动态规划最重要的还是找到递归式子,找准问题和子问题的关系。找到之后就能填表了,表填好了就可以回溯得到解了。