数据挖掘 发表于 2013-2-28 18:25:48

自动找出分子中环的编号

我2月8号来上海,今天28号,前前后后化了我将近20天的时间,才搞定了自动识别mol2中环的编号。
需要一些图论方面的知识,生物专业出身,数学和计算机底子薄。希望能借此认识更多的同行人,进行交流。

图论无向图中的闭合回路http://cn.qzs.qq.com/ac/b.gif                                           
1前提:对一连通分量P,将其用邻接矩阵表示法来表示
0代表不连通,1 代表连通
2 用广度优先算法求出连通分量P的支撑树(即生成树)
生成树:是一个极小连通子图,它含有图中全部顶点,但只有n-1条边。
    由深度优先搜索遍历得到的生成树,称为深度优先生成树。
由广度优先搜索遍历得到的生成树,称为广度优先生成树。
见下页无向图G7的两种生成树。

http://b121.photo.store.qq.com/psb?/V13dY5Qg2PAdhc/r2mh1a6WAO72NE*np.GH9mjSQ0BkkRXBMpeTBqKAFCc!/b/dA3JJkgwLgAA&bo=qgDQAAAAAAADAF8!http://b126.photo.store.qq.com/psb?/V13dY5Qg2PAdhc/rV.iwbJn3eM3qmDYPQvGLr23brPPuQsps6QlddY11HY!/b/dBBzKUuwBgAA&bo=gAHlAAAAAAADAEE!
这里我们使用的是广度优先算法,所以我们得到的是图(b).广度优先生成树BFS
通过BFS算法,把边的权重改为-1;

------------------------------------------------------------------------------------------------------------------------------------------------------------------------

结果当然要放到正文里了。
------------------------------------------------
下面是连通信息
1   2
1   5
1   8
2   3
2    13
3   6
3    14
4   6
4   7
4    15
5   7
5    11
7   9
8    16
9    10
9    18
10    12
10    17
11    12
11    21
12    22
13    16
13    23
13    24
16    25
16    26
17    20
17    27
18    19
18    28
19    20
19    29
20    30
------------------------------------------------------

这里是连通信息表达的分子式

http://b115.photo.store.qq.com/psb?/V13dY5Qg2PAdhc/NtwAVceO7LUjwsj*Y*MhbCn9kD0HFPKIqmz*tWnbsmQ!/b/dMi3jkStCgAA&bo=IAOxAQAAAAABALQ!

-----------------------------------------------------------------------------
这里是我用perl写的脚本寻找环的编号,

#!/usr/bin/perl -w
use strict;

#N=30E=3330个点, 33 条边

#构建邻接矩阵
my $x;
my $y;
my @matrix;
for($x=1;$x<=30;$x++)
{
for($y=1;$y<=30;$y++)
{
$matrix[$x][$y]=0;
}
}
open FH,"F:\\C\\bond";
while(<FH>)
{
/(\d+)\s+(\d+)/;
$matrix[$1][$2]=1;
$matrix[$2][$1]=1;
}
close FH;
#for($x=1;$x<=30;$x++)
#{
# for($y=1;$y<=30;$y++)
# {
# print "$matrix[$x][$y]";
# }
# print "\n";
#}
#打印邻接矩阵


my @color;      #对每个顶点u∈V,其色彩存储于变量color中.
my @parent;    #结点u的父母存于变量π=parent中
                #如果u没有父母(例如u=s或u还没有被检索到),则 π=NIL

##初始化
#把所有的点初始化为白色
for(1..30)
{
$color[$_]="w";
}
my @box;#用来存放灰色的顶点




$parent=0;   #表示1号原子没有父母,也就是说任命 1号原子为root,祖宗
$color="g";    #把1号原子的颜色变为灰色ggrey,
                #黑色bBlack;   灰色grey    g
               
push @box,"1";#把灰色的1号原子放到盒子中
   
   
               
#写一个子例程,用来返回一个原子的相连原子
#print join "\n",&adj(2);
subadj
{
my @adjatoms=();
my $adj_atom=$_;
open FH,"F:\\C\\bond";
while(<FH>)
{
if($_=~/\b$adj_atom\b\s+(\d+)/ || $_=~/(\d+)\s+\b$adj_atom\b/)
{
   push @adjatoms,$1;
}
}
return @adjatoms ;
}

my @adj;                #Adj表示图中和u相邻的节点集合
while(@box!=0)          #当盒子中没有灰色的原子时,停止运行
{

my $greyatom=$box;
my @atoms=&adj($greyatom);
# print "\@atoms[$greyatom] ";
# print join "->",@atoms;
# print "\n";
foreach my $atom(@atoms)
{
   if($color[$atom] eq "w")
   {
    $color[$atom]="g";
    $parent[$atom]=$greyatom;
    $matrix[$greyatom][$atom]=-1;
   # print "$greyatom $atom-1\n";
   
   
    push @box,$atom;
   }
}
$color[$greyatom]="b";

#print "$greyatom ";打印黑色的点也就是生生成树
shift @box;
#print join "--",@box;
#print "\n";
}


#2. 在邻接矩阵01 -1中寻找权值不是-1的边(当然也不是0,
#如果是无权图,就应该找值为1的边),
#假定该边连接的是节点i和j。将其边的权值改为-1;

my @dealnums;


for($x=1;$x<=30;$x++)
{
for($y=1;$y<=30;$y++)
{
if($matrix[$x][$y]==1)
   {
    $matrix[$x][$y]=-1;
    push @dealnums,[($x,$y)];
   
   
    print"这就我要用dfs处理的$x$y\n";
   
   
#    my $start;
#      my $end;
#      my @path;
#      $start=$x;                     
#      $end=$y;                        
#      @path=();                     
#      push @path,$start;            
#      
#      for(1..30)
#      {
#      $visit[$_]='n';#n代表没有被访问过,v代表被访问过了
#      }                  
#                                    
#      $visit[$start]='v';            
#                                    
#                                    
#      @pathsij=&dfs($start,$end,\@path);
#                                    
#                           
#      foreach my$ref(@pathsij)         
#      {                                 
#       print join"->", @{$ref};         
#       print "\n";                     
#      }                                 
                                       

   }
}

}
our @pathsij;
our @visit;
foreach my $ref(@dealnums)
{

my @numsxy=@{$ref};
my $start;
my $end;
my @path;
$start=$numsxy;
$end=$numsxy;
@path=();
push @path,$start;

for(1..30)
{
   $visit[$_]='n';#n代表没有被访问过,v代表被访问过了
}               
$visit[$start]='v';
&dfs($start,$end,\@path);
}
my %hash;
foreach my$ref(@pathsij)
{

if(2<@{$ref} && @{$ref}<=8)
{
my @nums=sort @{$ref};
my $key=join('.',@nums);
$hash{$key}=1;
#print join"->", sort @{$ref};
#print "\n";
}
}
print join "\n",my @keyss=keys(%hash);

############深度遍历的子例程dfs

#写一个子例程,用来返回一个原子的相连原子
#print join "\n",&adj(2);

subdfs
{
my $begin=$_;
my $terminal=$_;
my @dfsdots=@{$_};
#找出$begin的邻接顶点




if($begin==$terminal)
{
   #print "\n......";
   #print join "--", @dfsdots;
   #print "……………………\n";
   push @pathsij,\@dfsdots;
   #return@dfsdots ;
}
else
{

my @adjbeginatoms=&adj($begin);
#print "$begin@adjbeginatoms \n";

foreach my $atom(@adjbeginatoms)
{
#print "^^$atom^^\n";
   
   #寻找新的出发点 (Vi,Vj)∈E,且Vj未访问过,故Vj为新出发点
   if($visit[$atom] ne 'v')
   {
      push @dfsdots,$atom;
      $visit[$atom]='v';
            
      &dfs($atom,$terminal,\@dfsdots);
         # pop @dfsdots;
#push pop 同时放到内部和外部理论上是可行的
#内部的会减少操作,所以我选择了内部

$visit[$atom]='0';
   pop @dfsdots;   
         
   }
   
   }

}


}

----------------------------------------------------

这里是输出结果

10.11.12.5.7.9
1.2.3.4.5.6.7
10.17.18.19.20.9
1.13.16.2.8


----------------------------------------------------------------------------------------------------------------------------------
这里是我参考的伪代码

http://blog.csdn.net/lzrzhao/article/details/2175787

我也把这个伪代码贴出来。
2.对于每一个连通分量,单独计算其环的个数,
则无向图G的总环数即为各连通分量环数总和

前提:对一连通分量P,将其用邻接矩阵表示法来表示
(对于无权图可以用1表示)

1. 用广度优先算法求出P的支撑树(即生成树),在求支撑树的过程中,
用 -1表示被加入支撑树中的边。;
2. 在邻接矩阵01 -1中寻找权值不是-1的边(当然也不是0,
如果是无权图,就应该找值为1的边),
假定该边连接的是节点i和j。将其边的权值改为-1;
针对矩阵的
二重循环可以搞定

3. 采用深度优先遍历算法求出从顶点i到顶点j之间所有简单路径
这里应该也是针对图的。


(注意给每个顶点赋不同权值。
例如-1,0,1分别表示未遍历,
已经遍历但还有相邻结点未遍历完,
已经遍历而且相邻结点已遍历完。
这样做主要是为了防止回溯到上一已访问过的结点。);

3. 根据生成树的定义,在生成树上每增加一条边,就会有一个回路。
在生成树上寻找i和j的路径。将该简单路径与边(i,j)连接即得环。
输出该环;
4. 继续在邻接矩阵中寻找权值不是-1的边,
假定该边连接的两顶点是v和w。将其边的权值改为-1;
5. 求出从顶点i到顶点j之间的所有简单路径;
6.分别将所求出的简单路径与边(i,j)连接即得环,输出该环;
7.重复执行步骤4-7,直到在邻接矩阵中没有权值是-1的边为止。
-----------------------------------------------------------------------------------------------------------
伪代码的评价
只能迂回的解决问题,它是找所有的环,不排除环中环;
ij 找到的回路
和vw找到的回路
可能会有交叉,而且找到的回路不都是简单回路
7.重复执行步骤4-7,直到在邻接矩阵中没有权值是-1的边为止。

这里写错了吧,应该是没有1的边吧

---------------------------------------------------------------------

不管怎么样,非常感谢这位仁兄,把他的伪代码贴到网上,也感谢鱼儿老师,陪我解读了一些代码。
以及一直默默帮助我的gyh,huihui,corephi,当然师兄给我的动画也很形象。以及最近认识的pp.

感谢一直帮助我的perl qq群。


希望高手继续指点,提供更好的算法。





墨竹晓风 发表于 2013-3-1 08:53:35

厉害!
页: [1]
查看完整版本: 自动找出分子中环的编号