求教字符串快速比较问题

luffy.deng

回复 19# HyryStudio

回复  hulnglei
如果你希望在100万个和100万个之间任何一对进行比较的话，需要考虑一下算法。例如可以将其中的一组建立一棵树，字符串的每个字符是一个节点，它的子节点是下一个字符，因此树的深度是字符串的长度。然后在此树中对另外一组的每个字符串进行搜索，当发现有两个节点不同时就可以停止搜索。这时可以省略后续所有节点的搜索，应该可以成数量级的加快计算速度。

如果你的字符串中的字符只有ATCG四个字符的话，并且长度是31的话，可以考虑使用两个32bit的长整型数保存其状态。这样可以通过异或位操作快速比较，当然还需要一个数比特1个数的函数，这些在C语言级别做的话应该是非常快的。
HyryStudio 发表于 2010-06-02 20:52

这个很牛，不妨实现出来看一下?
关于字符串比较( 比如 c的strcmp)大体的实现都是这样

1.   
   
2. while (*s1 != 0 && *s1 == *s2)
   
3.     s1++, s2++;
   
4. 
   
5.   if (*s1 == 0 || *s2 == 0)
   
6.     return (unsigned char) *s1 - (unsigned char) *s2;
   
7.   return *s1 - *s2;
   
8.   

复制代码

内核里的实现原理跟上面一样，只是用了内联汇编

1. extern inline int strcmp(const char * cs,const char * ct)
   
2. {
   
3. register int __res __asm__("ax");
   
4. __asm__("cld\n"
   
5.         "1:\tlodsb\n\t"
   
6.         "scasb\n\t"
   
7.         "jne 2f\n\t"
   
8.         "testb %%al,%%al\n\t"
   
9.         "jne 1b\n\t"
   
10.         "xorl %%eax,%%eax\n\t"
    
11.         "jmp 3f\n"
    
12.         "2:\tmovl $1,%%eax\n\t"
    
13.         "jb 3f\n\t"
    
14.         "negl %%eax\n"
    
15.         "3:"
    
16.         :"=a" (__res):"D" (cs),"S" (ct):"si","di");
    
17. return __res;
    
18. }

复制代码

hulnglei

回复  hulnglei
如果你希望在100万个和100万个之间任何一对进行比较的话，需要考虑一下算法。例如可以将其 ...
HyryStudio 发表于 2010-06-02 20:52

   谢谢，我当时听C语言数据结构的没好好听。C语言忘的一干二净。这树怎么做我真还得好好学学。现在项目赶，暂时是没工夫了。

hulnglei

回复  HyryStudio


这个很牛，不妨实现出来看一下?
关于字符串比较( 比如 c的strcmp)大体的实现都是这 ...
luffy.deng 发表于 2010-06-02 22:18

呵呵！后来发现自己还是搞错了。不用树来算的话。还是很慢。算还是得一个多月。这肯定是不行的了。我只能找找有没有用树算的软件。树的话只能请教我朋友写下了。谢谢各位！

HyryStudio

回复 23# hulnglei

用树写其实也不难，不过最好还是用C语言写，才能够达到飞快，不过这个比较费时间了。我想了一个很好的办法，你看是否可行：

假设A,B两个组，每个180万个长度为31的每位有ATCG等四种可能的字符串。
把B组做成一个集合set，或者字典，关键字是字符串，值是它们的名字。

对于A组中的每一个字符串，在B中找是否有和它相差一个元素的。
由于字典只能找完全一样的，因此具体这样做：

因为你只找相差一个字符的，因此对于字符串中的每个字符，都修改为其它三个，这样31*3 = 93，也就是每个字符串搜索93次。也就是穷举和此字符串相差一个的所有可能情况，然后都到字典中去找：
例如对于

AC....
搜索
TC...., CC..., GC..., AA..., AT..., AG...,


这样搜索空间从180万*180万降低到180万*93，应该能在有限的时间之内完成。

另外，先用2个32位的整数对数据进行编码再做成字典的话，能够节省不少内存。不过这一步不做也可，只要你内存足够大。

Kabie

楼上这解法很棒啊。。。

hulnglei

回复  hulnglei

用树写其实也不难，不过最好还是用C语言写，才能够达到飞快，不过这个比较费时间了。我 ...
HyryStudio 发表于 2010-06-03 19:41

   我试完。几分钟就算完了！太爽了。终于真的可以拿来用了！太感谢了！

luffy.deng

回复  hulnglei

因为你只找相差一个字符的，因此对于字符串中的每个字符，都修改为其它三个，这样31*3 = 93，也就是每个字符串搜索93次。也就是穷举和此字符串相差一个的所有可能情况，然后都到字典中去找：
例如对于

AC....
搜索
TC...., CC..., GC..., AA..., AT..., AG...,


这样搜索空间从180万*180万降低到180万*93，应该能在有限的时间之内完成。 ...
HyryStudio 发表于 2010-06-03 19:41

这个办法很很牛X,我简单试了一下在我的机器上只跑了14分钟多。

1. s ='AATTCTAACATGAAAGTAGGAAAGATGTCAC'
   
2. d={'AATTCTAACATGAAAGTAGGAAAGATGTCAC':0}
   
3. def genstr(s,d):
   
4.     n=31
   
5.     while n>0:
   
6.         for i in 'ATCG':
   
7.             if i !=s[n-1:n]:
   
8.                 d.has_key(s[0:n-1]+i+s[n:])
   
9.         n-=1
   
10. t=1800000
    
11. while t>0:
    
12.     print t
    
13.     genstr(s,d)
    
14.     t-=1

复制代码

GOD_Father

思想无极限。。受益良深

HyryStudio

回复 24# HyryStudio

用Python写了一个用搜索树实现的，结果发现当树很平衡很大时，每个字符串的搜索量基本上和93差不了多少。穷举法对于Python来说还是最快的，因为dict是C语言级别实现的。

hulnglei

回复  hulnglei

用树写其实也不难，不过最好还是用C语言写，才能够达到飞快，不过这个比较费时间了。我 ...
HyryStudio 发表于 2010-06-03 19:41

   这样我顺便再问下，有兴趣的高手可以帮忙想下。如果两个字符串不是只有一个字符不同。有两个，三个，或者四个呢？ 在生物上面单个字符可以解释为SNP， 2-4个也是可以解释为单倍型。也是有很用。

1.                 fo = file("out_3snp.txt",'w')
   
2.                 for keya,vala in self.adict.items():
   
3.                         if vala in self.bdict:
   
4.                                 continue
   
5.                         else:
   
6.                                 for i1 in range(5, 31):
   
7.                                         for j1 in 'ACGT':
   
8.                                                 if vala[i1] == j1:
   
9.                                                         continue
   
10.                                                 else:
    
11.                                                         b1= vala[0:i1]+j1+vala[i1+1:31]
    
12. 
    
13.                                                 for i2 in range(5,31):
    
14.                                                         if i2==i1:
    
15.                                                                 if (b1 not in self.cdict) and (b1 in self.bdict):
    
16.                                                                         print >> fo, keya, self.bdict[b1][0], vala[:31], i1+1, vala[i1+32]
    
17.                                                                 break
    
18.                                                         else:
    
19.                                                                 for j2 in 'ATCG':
    
20.                                                                         if vala[i2] == j2:
    
21.                                                                                 continue
    
22.                                                                         else:
    
23.                                                                                 b2 = b1[0:i2]+j2+b1[i2+1:31]
    
24. 
    
25.                                                                         for i3 in range(5, 31):
    
26.                                                                                 if i3==i1 or i3==i2:
    
27.                                                                                         if (b2 not in self.cdict) and (b2 in self.bdict):
    
28.                                                                                                 print >> fo, keya, self.bdict[b2][0], vala[:31], i1+1, vala[i1], j1, vala[i1+32],
    
29. i2+1, vala[i2], j2, vala[i2+32]
    
30.                                                                                         break
    
31.                                                                                 else:
    
32.                                                                                         for j3 in 'ACGT':
    
33.                                                                                                 if vala[i3] == j3:
    
34.                                                                                                         continue
    
35.                                                                                                 else:
    
36.                                                                                                         b3 = b2[0:i3]+j3+b2[i3+1:31]
    
37. 
    
38.                                                                                                 if (b3 not in self.cdict) and (b3 in self.bdict):
    
39.                                                                                                         print >> fo, keya, self.bdict[b3][0], vala[:31], i1+1, vala[i1], j1, vala[i1+32], i2+1, vala[i2], j2, vala[i2+32], i3+1, vala[i3], j3, vala[i3+32]
    
40. -- INSERT --                                                                              

复制代码