生物分子模拟论坛

 找回密码
 我想注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

查看: 1294|回复: 1

[Perl] 自动找出分子中环的编号

[复制链接]
发表于 2013-2-28 18:25:48 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,下载更多分子模拟资源。

您需要 登录 才可以下载或查看,没有帐号?我想注册

x
我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群。


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





发表于 2013-3-1 08:53:35 | 显示全部楼层
厉害!
您需要登录后才可以回帖 登录 | 我想注册

本版积分规则

QQ|分迪科技|小黑屋|手机版|Archiver|生物分子模拟论坛 ( 蜀ICP备14009200号-3

GMT+8, 2018-8-21 15:46 , Processed in 0.254884 second(s), 29 queries .

Powered by Discuz! X3.3

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表