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

作者:bibodeng 发布于:2012-4-24 8:08 Tuesday 分类:编程技术

【前言】
对于一个问题,我们如果可以枚举所有的解,那么这个解空间我们是知道的。那么如何在解空间里面找到最优解呢?这时有一个非常好的方法,从底向上地构造整个解,而每一步都是从地层寻求最优解,这样就能保证在最终得到的一定是最优解。这就是最优子结构,有这种结构的问题,通常都可以用动态规划的办法来寻求最优解。而且它是从小规模(子问题)到大规模问题的构造,而常常这样的解法能够用一张表直观地表现出来。表中的元素是一个表达式的某个特定值,这个表达式表现的是问题和子问题的关系,也就是如何在子问题中的解中寻找最优的关系,这样的关系在例子中会非常地明了。
 
【装配流水线】
 
往往最经典的算法书里面都会讲最经典的“装配流水线”问题。因为它相对来说比较简单,一个问题,只有两个个子问题,那就是两条流水线选哪条路径。即使是最简单的例子也会有一大通的表达式搞得头有点晕,不过多看几遍应该是可以克服的。
假设一个工厂有两条装配线:装配线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] 比较特殊,因为一开始只能从一条路径进入,用时为e(i 为0或1)

Fi[0] = e+ 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

标签: 算法 编程 DNA 动态规划

发表评论:

Powered by emlog 京ICP备16017775