自动找出分子中环的编号
我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群。
希望高手继续指点，提供更好的算法。
厉害！
页:
[1]