做了一下中科大预赛C题，题目在http://202.38.68.79/~mcm/problem.pdf
现把程序贴在下面：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % %
% MAIN %
% % % %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%最短DNA子序列(中科大MCM预赛题第三题)
%此主函数提供两个功能，1.给出你已获得的DNA片段，合成最短的DNA单链;
%2.你给出所要测试的DNA片段的元素数及片段个数，程序随机生成相应片段，再合成；
%simula.m就是完成随机生成功能模块的。
%str_eva.m功能是评价函数值，即计算总共有多少匹配。
%dou_cat.m的功能是计算两个片段有多少匹配。
clear
clc
seg='';
disp('Please input segments you have.And you should use the format of')
disp('"seg=zeros(n,m)" first to create the argument "seg" which memory all the segments')
disp('in which the argument n and m part represent the number of the segments ')
disp('and the number of columns in every segment.Secondly,one segment will be')
disp('made up to one row in argument seg.Then please press "Return" to continue.')
disp('And if you want to create the segments stochasticly.Please press "Return" to continue.')
keyboard
if isempty(seg)
disp('Please input the argument n which represent the number of ')
disp('the segments.(and you should use the format of "n=?")')
disp('Please input the argument m which represent the number of ')
disp('the columns in a segment.(and you should use the format of "m=?")')
disp('Then you should press "Return" to continue.')
keyboard
end
if ~isempty(seg)
seg=char(seg);
[DNA,len_match,len_seg]=sho_dna(seg)
else
if (isempty(n)|isempty(m))
disp('n and m must not be empty!')
else
seg=simula(n,m);
seg=char(seg);
[DNA,len_match,len_seg]=sho_dna(seg)
end
end
function [DNA,len_mat,len_seg]=sho_dna(seg)
%用模拟退火算法实现
%order指各个串的次序，seg保存所有片段
%初始化变量
order_m=[1:size(seg,1)];
order_pro=order_m;
order_save=order_m; % 意义同len_m,len_pro,len_save
len_m=0;len_pro=0;len_save=0; %len_m某一次的,len_pro过程中的，len_save保存最大值的
seg_m='';seg_pro='';seg_save='';
t0=5;tf=3;a=0.83;t=t0;t_mid=25;
while t>=tf
for i=1:size(seg,1)
%产生随机扰动
%若温度高于某个值，则产生不止一次扰动
if t>t_mid
m=2;
else
m=1;
end
for j=1:m
k=ceil((size(seg,1)-1).*rand)+1;
temp=order_m(k);
order_m(k)=order_m(k-1);
order_m(k-1)=temp;
end
%评价
[seg_m,len_m]=str_eva(seg,order_m);
%判断是否接受
if len_pro<len_m
len_pro=len_m;
order_pro=order_m;
seg_pro=seg_m;
else
rand_jud=rand;
if rand_jud>=exp(-(len_pro-len_m)./t)
len_pro=len_m;
order_pro=order_m;
seg_pro=seg_m;
else
order_m=order_pro;
end
end
%保存最大的
if len_m>len_save
len_save=len_m;
order_save=order_m;
seg_save=seg_m;
end
end
t=t.*a;
end
DNA=seg_save;
len_mat=len_save;
len_seg=size(seg_save,2);
function [str_res,count]=dou_cat(str1,str2)
%str1要求不短于str2
n=size(str1,2);
m=size(str2,2);
str22=str2(m:-1:1);
q=0;
for k=m:-1:1
if sum(str1(n-k+1:n)==str2(1:k))==k %还差那个p编号的问题
if m==k
stra='';
strb='';
strm=str2(1:m);
count=m;
q=1;
break
else
stra=str1(1:n-k);
strb=str2(k+1:m);
strm=str2(1:k);
count=k;
q=1;
break
end
elseif sum(str2(m-k+1:m)==str1(1:k))==k
if k==m
stra='';
strb='';
strm=str2(1:m);
else
stra=str2(1:m-k);
strb=str1(k+1:n);
strm=str1(1:k);
count=k;
q=1;
break
end
elseif sum(str1(n-k+1:n)==str22(1:k))==k
if m==k
stra='';
strb='';
strm=str1(1:m);
count=m;
q=1;
break
else
stra=str1(1:n-k);
strb=str22(k+1:m);
strm=str22(1:k);
count=k;
q=1;
break
end
elseif sum(str22(m-k+1:m)==str1(1:k))==k
stra=str22(1:m-k);
strb=str1(k+1:n);
strm=str1(1:k);
count=k;
q=1;
break
end
end
if q==1
str_res=strcat(stra,strm);
str_res=strcat(str_res,strb);
else
str_res=strcat(str1,str2);
count=0;
end
function [str_gen,count_gen]=str_eva(seg,order)
%评价函数值
seg_mo=seg(order(1), ;
count=0;
count_gen=0;
for i=2:size(seg,1)
[seg_mo,count]=dou_cat(seg_mo,seg(order(i),);
count_gen=count_gen+count;
end
str_gen=seg_mo;
function seg=simula(m,n) %随机产生m行,n列的片段
seg=zeros(m,n);
for i=1:m
for j=1:n
temp=ceil(rand.*4);
switch temp
case 1
seg(i,j)='A';
case 2
seg(i,j)='C';
case 3
seg(i,j)='G';
case 4
seg(i,j)='T';
end
end
end
[此贴子已经被作者于2004-1-30 18:38:07编辑过]
|
发表于 2004-1-24 07:49:48
楼主
发表于 2004-1-31 08:08:08