【SCMS】最短公共母序列的代码实现

大家好，欢迎来到IT知识分享网。

前言

其实说起初衷，是之前吾辈的学校开展了个数模比赛，我去帮别人当帮手，当时的B题是基因拼接，D题是污染程度评定

作为ACMer，我当然喜欢B题这种像ACM竞赛题一样感觉的东西，但无奈人家才是参(da)赛(ye)者(a)，他们选了D，嘛，其实据说那次我做的也还是不错。

回到正题，于是我就开始考虑这种算法改如何实现呢，反正现在是学生，在CNKI上看论文不要钱，到处看看咯。

看到个论文，用dp实现生物序列最短公共超序列，感觉靠谱，要不然就试试呗，拿来当模板用，说不定哪天就出道题呢~

给他起了个名儿

SCMS = Shortest Common Mother-Sequence (最短公共母序列)

建模

可以将题意理解为ACM格式，则为如下情况：

数据组数为T，【给出T】

每组数据中有n条基因序列，【给出N】

以下则为n行，每行以一个字符串表示一条基因序列的组成【给出str_i】

【每条基因序列均为母序列中的一个子片段，此处的子片段指的是子串而非子序列】（给自己降低一点难度嘛不要在意）

要求输出满足所有这些基因序列的母序列。

代码实现(C/C++)

#include<stdio.h> #include<string.h> struct str { char letter[1600]; int len; }ans,a[20]; struct { struct str temp; int x,y; }add[500]; int f[1<<12][15],father[1<<12][15],reduce[20][20],cost[20][20],n,m,max,Test,s; struct str ADD(struct str ans,int k,int reduce) { int i; for(i=ans.len;i<ans.len+a[k].len-reduce;i++) ans.letter[i]=a[k].letter[reduce+i-ans.len]; ans.letter[i]=0; ans.len=ans.len+a[k].len-reduce; return ans; } void quick(int l,int r) { int i=l,j=r; struct str x=add[(i+j)>>1].temp; while(i<=j) { while(strcmp(add[i].temp.letter,x.letter)<0) i++; while(strcmp(add[j].temp.letter,x.letter)>0) j--; if(i<=j) { add[0]=add[i],add[i]=add[j],add[j]=add[0]; i++,j--; } } if(i<r) quick(i,r); if(j>l) quick(l,j); } void READY() { int i,j,k,r; for(i=1;i<=n;i++) for(j=1;j<=n;j++) for(k=0;k<a[i].len;k++) { for(r=0;k+r<a[i].len&&a[i].letter[k+r]==a[j].letter[r];r++); if(k+r==a[i].len) { reduce[j][i]=r; break; } } s=0; for(i=1;i<=n;i++) for(j=1;j<=n;j++) { s++; add[s].temp=ADD(a[i],j,reduce[j][i]),add[s].x=i,add[s].y=j; } quick(1,s); for(i=1;i<=s;i++) cost[add[i].x][add[i].y]=i; } void read() { int i,j,k,r; scanf("%d",&n); memset(reduce,0,sizeof(reduce)); for(i=1;i<=n;i++) { scanf("%s",a[i].letter); a[i].len=strlen(a[i].letter); } for(i=1;i<=n;i++) for(j=1;j<=n;j++) { if(i==j) continue; for(k=0;k<a[j].len;k++) { for(r=0;r+k<a[j].len&&r<a[i].len&&a[j].letter[k+r]==a[i].letter[r];r++); if(r==a[i].len) break; } if(k<a[j].len) a[i].letter[0]=0; } for(i=n;i>=1;i--) if(a[i].letter[0]==0) { for(j=i;j<n;j++) a[j]=a[j+1]; n--; } m=(1<<n)-1; } void DP() { int i,j,k; for(i=0;i<=m;i++) for(j=1;j<=n;j++) f[i][j]=father[i][j]=0; for(i=1;i<=m;i++) for(j=1;j<=n;j++) { int add=1<<(j-1); if((i&add)==0) for(k=1;k<=n;k++) if((i&(1<<(k-1)))!=0) { if(father[i+add][j]==0) { f[i+add][j]=f[i][k]+reduce[k][j],father[i+add][j]=k; continue; } if(f[i+add][j]<f[i][k]+reduce[k][j]) f[i+add][j]=f[i][k]+reduce[k][j],father[i+add][j]=k; else if(f[i+add][j]==f[i][k]+reduce[k][j]&&cost[j][father[i+add][j]]>cost[j][k]) father[i+add][j]=k; } } } void solve() { max=-1; int i,j,now; for(i=1;i<=n;i++) { if(max<f[m][i]) { max=f[m][i],now=m; ans=a[i]; for(j=i;father[now][j]!=0;) { ans=ADD(ans,father[now][j],reduce[father[now][j]][j]); int t=j; j=father[now][j],now-=(1<<(t-1)); } } else if(max==f[m][i]) { struct str temp=a[i]; now=m; for(j=i;father[now][j]!=0;) { temp=ADD(temp,father[now][j],reduce[father[now][j]][j]); int t=j; j=father[now][j],now-=(1<<(t-1)); } if(strcmp(ans.letter,temp.letter)>0) ans=temp; } } } int main() { int T=0; scanf("%d",&T); for(int i=1;i<=T;i++) { read(); READY(); DP(); solve(); printf("%s\n",ans.letter); } return 0; } 

Reference

参阅的论文如下，关于证明过程和各种条件约束详见论文：

 《两个生物序列最短公共超序列的动态规划算法》

黄永莲1, 孙世军2

 (1.湛江师范学院生物科学与化学学院,广东湛江; 

摘要:将动态规划算法用于两个生物序列的最短公共超序列的计算。

计算过程分为两个算法,第一个算法计算两个序列的所有前缀的最短公共超序列的长度,

并存放在一个矩阵中,第二个算法利用前面所得的矩阵,找到两个序列的最短公共超序列。

关键词:生物序列;动态规划算法;超序列;最短公共超序列