自动找出分子中环的编号

数据挖掘

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

图论  无向图中的闭合回路

                               
登录/注册后可看大图

                                            
1前提：对一连通分量P，将其用邻接矩阵表示法来表示
0代表不连通，1 代表连通
2 用广度优先算法求出连通分量P的支撑树（即生成树）
生成树：是一个极小连通子图，它含有图中全部顶点，但只有n-1条边。
    由深度优先搜索遍历得到的生成树，称为深度优先生成树。
由广度优先搜索遍历得到的生成树，称为广度优先生成树。
见下页无向图G7的两种生成树。

                               
登录/注册后可看大图

                               
登录/注册后可看大图

这里我们使用的是广度优先算法，所以我们得到的是图（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
------------------------------------------------------

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

                               
登录/注册后可看大图

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

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

#N=30  E=33  30个点， 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[1]=0;   #表示1号原子没有父母，也就是说任命 1号原子为root，祖宗
$color[1]="g";    #把1号原子的颜色变为灰色g  grey，
                #黑色b  Black;   灰色grey    g
               
push @box,"1";  #把灰色的1号原子放到盒子中
   
   
               
#写一个子例程，用来返回  一个原子的相连原子
#print join "\n",&adj(2);
sub  adj
{
my @adjatoms=();
my $adj_atom=$_[0];
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[0];
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[0];
  $end=$numsxy[1];
  @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);

sub  dfs
{
my $begin=$_[0];
my $terminal=$_[1];
my @dfsdots=@{$_[2]};
#找出$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群。


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

墨竹晓风

厉害！