在相同长度的两个字符串之间找到不匹配位置的最快方法

Question

我有数百万对相同长度的字符串，我想比较并找到它们 不匹配的位置。

例如，对于每个$str1 和$str2，我们要查找不匹配 位置与$str_source:

$str_source = "ATTCCGGG";

$str1       = "ATTGCGGG"; # 1 mismatch with Str1 at position 3 (0-based)
$str2       = "ATACCGGC"; # 2 mismatches with source at position  2 and 7

有没有快速的方法。目前我有循环的 C 风格方法 使用“substr”函数在两个字符串中的每个位置。但是这种方法非常慢。

my @mism_pos;
for $i (0 .. length($str_source)) {
   $source_base = substr($str_source,$i,1);
   $str_base    = substr($str2,$i,$1);

  if ($source_base ne $str_base) {
     push @mism_pos,$i;
  }

}

Answer 1

Inline::C

计算很简单，用Inline::C （阅读 perldoc Inline::C-Cookbook 和 perldoc Inline::C 获取文档）：

use Inline C => << '...';                                                       
  void find_diffs(char* x, char* y) {                                           
    int i;                                                                      
    Inline_Stack_Vars;                                                          
    Inline_Stack_Reset;                                                         
    for(i=0; x[i] && y[i]; ++i) {                                               
      if(x[i] != y[i]) {                                                        
        Inline_Stack_Push(sv_2mortal(newSViv(i)));                              
      }                                                                         
    }                                                                           
    Inline_Stack_Done;                                                          
  }                                                                             
...                                                                             

@diffs= find_diffs("ATTCCGGG","ATTGCGGG");  print "@diffs\n";                   
@diffs= find_diffs("ATTCCGGG","ATACCGGC");  print "@diffs\n";                   

这是该脚本的输出：

> script.pl 
3
2 7

PDL

如果您想在 Perl 中快速处理大量数据，请学习 PDL (Documentation)：

use PDL; 
use PDL::Char;                                                                  
$PDL::SHARE=$PDL::SHARE; # keep stray warning quiet 

my $source=PDL::Char->new("ATTCCGGG");                                          
for my $str ( "ATTGCGGG", "ATACCGGC") {                                         
  my $match =PDL::Char->new($str);                                              
  my @diff=which($match!=$source)->list;                                        
  print "@diff\n";                                                              
}

（与第一个脚本的输出相同。）

注意事项：我在基因组数据处理中使用 PDL 非常愉快。连同对存储在磁盘上的数据的内存映射访问，可以快速处理大量数据：所有处理都在高度优化的 C 循环中完成。此外，对于PDL 中缺少的任何功能，您可以通过Inline::C 轻松访问相同的数据。

但是请注意，创建一个 PDL 向量非常慢（恒定时间，对于大型数据结构是可以接受的）。因此，您宁愿一次性创建一个包含所有输入数据的大型 PDL 对象，也不愿遍历单个数据元素。

Answer 2

那些看起来像基因序列。如果字符串都是 8 个字符，并且可能的代码域是 (A, C, G, T)，您可能会考虑在处理数据之前以某种方式转换数据。那只会给你 65536 个可能的字符串，所以你可以专门化你的实现。

例如，您编写了一个方法，该方法接受一个 8 个字符的字符串并将其映射到一个整数。 Memoize 这样操作会很快。接下来，编写一个比较函数，给定两个整数，告诉你它们有何不同。在调用比较之前，您可以在一个合适的循环结构中调用它，并使用像 unless ( $a != $b ) 这样的数字相等测试 - 如果您愿意，可以将相同的代码短路。

Answer 3

听起来这可能是您的应用程序的性能关键部分。在这种情况下，您可能需要考虑编写一个 C 扩展方法来进行比较。

Perl 提供了XS 扩展机制，这使得这变得相当简单。

Answer 4

这是一个基准测试脚本，用于确定各种方法的速度差异。请记住，在第一次调用使用 Inline::C 的脚本时，在调用 C 编译器等时会有延迟。因此，运行一次脚本，然后进行基准测试。

#!/usr/bin/perl

use strict;
use warnings;

use Benchmark qw( cmpthese );

my ($copies) = @ARGV;
$copies ||= 1;

my $x = 'ATTCCGGG' x $copies;
my $y = 'ATTGCGGG' x $copies;
my $z = 'ATACCGGC' x $copies;

sub wrapper { 
    my ($func, @args) = @_;
    for my $s (@args) {
        my $differences = $func->($x, $s);
        # just trying to ensure results are not discarded
        if ( @$differences == 0 ) { 
            print "There is no difference\n";
        }
    }
    return;
}

cmpthese -5, {
    explode  => sub { wrapper(\&where_do_they_differ, $y, $z) },
    mism_pos => sub { wrapper(\&mism_pos, $y, $z) },
    inline_c => sub {
        wrapper(\&i_dont_know_how_to_do_stuff_with_inline_c, $y, $z) },
};

sub where_do_they_differ {
    my ($str1, $str2) = @_;
    my @str1 = split //, $str1;
    my @str2 = split //, $str2;
    [ map {$str1[$_] eq $str2[$_] ? () : $_} 0 .. length($str1) - 1 ];
}

sub mism_pos {
    my ($str1, $str2) = @_;
    my @mism_pos;

    for my $i (0 .. length($str1) - 1) {
        if (substr($str1, $i, 1) ne substr($str2, $i, 1) ) {
            push @mism_pos, $i;
        }
    }
    return \@mism_pos;
}

sub i_dont_know_how_to_do_stuff_with_inline_c {
    [ find_diffs(@_) ];
}

use Inline C => << 'EOC';
    void find_diffs(char* x, char* y) {
        int i;
        Inline_Stack_Vars;
        Inline_Stack_Reset;
        for(i=0; x[i] && y[i]; ++i) {
            if(x[i] != y[i]) {
                Inline_Stack_Push(sv_2mortal(newSViv(i)));
            }
        }
        Inline_Stack_Done;
    }
EOC

$copies = 1 的结果（在带有 AS Perl 5.10.1 的 Windows XP 上使用 VC++9）：

速率爆炸 mism_pos inline_c 爆炸 15475/s -- -64% -84% mism_pos 43196/s 179% -- -56% inline_c 98378/s 536% 128% --

$copies = 100 的结果：

速率爆炸 mism_pos inline_c 爆炸 160/s -- -86% -99% mism_pos 1106/s 593% -- -90% inline_c 10808/s 6667% 877% --

Answer 5

对于每个字符比较，您都对 substr 进行了 2 次调用，这可能会拖慢您的速度。

我会做一些优化

@source = split //,$str_source  #split first rather than substr
@base = split //, $str_base

for $i (0 .. length($str_source)) {
   $mism_pos{$1} = 1 if ($source[$i] ne $base); #hashing is faster than array push
}

return keys $mism_pos

Answer 6

比较字符串以找出差异的最快方法是将它们中的每个字节XOR，然后测试为零。如果我必须这样做，我只会用 C 编写一个程序来完成不同的工作，而不是编写 Perl 的 C 扩展，然后我会将我的 C 程序作为 Perl 的子进程运行。确切的算法将取决于字符串的长度和数据量。然而，这不会超过 100 行 C。事实上，如果你想最大限度地提高速度，可以用汇编语言编写一个对固定长度字符串的字节进行 XOR 并测试为零的程序。

Answer 7

一些经典的字符串比较优化：

最佳不匹配 - 从搜索字符串的 END 开始比较。例如如果您在开始时进行比较，则在 ABDABEABF 中搜索 ABC，您将一次沿着模式移动一个字符。如果你从头开始搜索，你将能够跳过三个字符

错误字符启发式 - 选择最不常见的字符并首先搜索。例如在英语中，“z”字符很少见，良好的字符串搜索功能将搜索“迷宫”并从第 3 个字符开始比较

Answer 8

我不知道它的效率如何，但是您总是可以对要匹配的两个字符串进行异或运算，并找到第一个不匹配的索引。

#! /usr/bin/env perl
use strict;
use warnings;
use 5.10.1;

my $str_source = "ATTCCGGG";

my $str1       = "ATTGCGGG";
my $str2       = "ATACCGGC";
my $str3       = "GTTCCGGG";

# this returns the index of all of the mismatches (zero based)
# it returns an empty list if the two strings match.
sub diff_index{
  my($a,$b) = @_;
  my $cmp = $a^$b;

  my @cmp;
  while( $cmp =~ /[^\0]/g ){ # match non-zero byte
    push @cmp, pos($cmp) - 1;
  }
  return @cmp;
}

for my $str ( $str_source, $str1, $str2, $str3 ){
  say '# "', $str, '"';
  my @ret = diff_index $str_source, $str;
  if( @ret ){
    say '[ ', join( ', ', @ret), ' ]';
  }else{
    say '#   match';
  }
}

# "ATTCCGGG"
#   match
# "ATTGCGGG"
[ 3 ]
# "ATACCGGC"
[ 2, 7 ]
# "GTTCCGGG"
[ 0 ]

---

通过B::Concise 运行它表明 CPU 昂贵的操作是作为单个操作码发生的。这意味着这些操作是在 C 中运行的。

perl -MO=Concise,-exec,-compact,-src,diff_index test.pl |
perl -pE's/^[^#].*? \K([^\s]+)$/# $1/' # To fix highlighting bugs

main::diff_index:
# 15:   my($a,$b) = @_;
1  <;> nextstate(main 53 test.pl:15) # v:%,*,&,$
2  <0> pushmark # s
3  <$> gv(*_) # s
4  <1> rv2av[t3] # lK/3
5  <0> pushmark # sRM*/128
6  <0> padsv[$a:53,58] # lRM*/LVINTRO
7  <0> padsv[$b:53,58] # lRM*/LVINTRO
8  <2> aassign[t4] # vKS
# 16:   my $cmp = $a^$b;
9  <;> nextstate(main 54 test.pl:16) # v:%,*,&,$
a  <0> padsv[$a:53,58] # s
b  <0> padsv[$b:53,58] # s
c  <2> bit_xor[t6] # sK                     <-----  Single OP -----
d  <0> padsv[$cmp:54,58] # sRM*/LVINTRO
e  <2> sassign # vKS/2
# 18:   my @cmp;
f  <;> nextstate(main 55 test.pl:18) # v:%,*,&,{,$
g  <0> padav[@cmp:55,58] # vM/LVINTRO
# 20:   while( $cmp =~ /[^\0]/g ){ # match non-zero byte
h  <;> nextstate(main 57 test.pl:20) # v:%,*,&,{,$
i  <{> enterloop(next->r last->v redo->j) # v
s  <0> padsv[$cmp:54,58] # s
t  </> match(/"[^\\0]"/) # sKS/RTIME        <-----  Single OP -----
u  <|> and(other->j) # vK/1
# 21:     push @cmp, pos($cmp) - 1;
j      <;> nextstate(main 56 test.pl:21) # v:%,*,&,$
k      <0> pushmark # s
l      <0> padav[@cmp:55,58] # lRM
m      <0> padsv[$cmp:54,58] # sRM
n      <1> pos[t8] # sK/1
o      <$> const(IV 1) # s
p      <2> subtract[t9] # sK/2
q      <@> push[t10] # vK/2
r      <0> unstack # v
           goto # s
v  <2> leaveloop # vK/2
# 24:   return @cmp;
w  <;> nextstate(main 58 test.pl:24) # v:%,*,&,{,$
x  <0> pushmark # s
y  <0> padav[@cmp:55,58] 
z  <@> return # K
10 <1> leavesub[1 ref] # K/REFC,1

Answer 9

我也想说“用 C 写”。

在那里，您可以使用优化，例如一次比较 4 个字符（作为 32 位整数）。

或者更改您的表示（4 个字母，对吗？）以使用 2 位表示基 (?)，这样您就可以一次比较 16 个字符。