Beefy Boxes and Bandwidth Generously Provided by pair Networks
Your skill will accomplish
what the force of many cannot
 
PerlMonks  

Hamming Distance Between 2 Strings - Fast(est) Way?

by monkfan (Curate)
on Oct 14, 2005 at 14:30 UTC ( #500235=perlquestion: print w/ replies, xml ) Need Help??
monkfan has asked for the wisdom of the Perl Monks concerning the following question:

My most revered monks,

The subroutine below compute the number of mismatches between two strings (they are always equal) - usually called Hamming Distance.

However my code below is painfully slow. I am terribly in need of a super fast way to do this. In particular it has to process millions of string pairs.

Thus I turn to you my brother monks, for illumination in this matter.
#!/usr/bin/perl -w use strict; my $s1 = 'AAAAA'; my $s2 = 'ATCAA'; my $s3 = 'AAAAA'; hd($s1,s2) # will give value 2 hd($s1,s3) # will give value 0 sub hd { #String length is assumed to be equal my ($k,$l) = @_; my $len = length ($k); my $num_mismatch = 0; for (my $i=0; $i<$len; $i++) { ++$num_mismatch if substr($k, $i, 1) ne substr($l, $i, 1); } return $num_mismatch; }
Update: Benchmark. Thanks so much everybody.
Rate Mine RJ BWK inman Mine 176987/s -- -65% -78% -80% RJ 504123/s 185% -- -38% -44% BWK 817943/s 362% 62% -- -9% inman 901871/s 410% 79% 10% -- *Note: Mine has already included the modification suggested by 'wazoox'.

Regards,
Edward

Comment on Hamming Distance Between 2 Strings - Fast(est) Way?
Select or Download Code
Re: Hamming Distance Between 2 Strings - Fast(est) Way?
by NateTut (Deacon) on Oct 14, 2005 at 14:39 UTC
      Dear NateTut,
      This module doesn't compute the exact Hamming Distance as I required.
      Instead it returns edit distance. Here is the example:
      perl -MText::LevenshteinXS -e ' $s1 = 'GGAAG'; $s2 = 'GAAGA'; $diff = distance($s1,$s2); print "$diff\n";'
      It prints: 2 instead of 3.

      Regards,
      Edward
Re: Hamming Distance Between 2 Strings - Fast(est) Way?
by BrowserUk (Pope) on Oct 14, 2005 at 14:51 UTC

    #!/usr/bin/perl -slw use strict; my $s1 = 'AAAAA'; my $s2 = 'ATCAA'; my $s3 = 'AAAAA'; print "$s1:$s2 hd:", hd( $s1, $s2 ); # will give value 2 print "$s1:$s3 hd:", hd( $s1, $s3 ); # will give value 0 sub hd{ length( $_[ 0 ] ) - ( ( $_[ 0 ] ^ $_[ 1 ] ) =~ tr[\0][\0] ) } __END__ P:\test>500235 AAAAA:ATCAA hd:2 AAAAA:AAAAA hd:0

    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
    "Science is about questioning the status quo. Questioning authority".
    The "good enough" maybe good enough for the now, and perfection maybe unobtainable, but that should not preclude us from striving for perfection, when time, circumstance or desire allow.
      Depending on the length and number of your nucleic acids, and the number of comparisons you are making it might be better to store and compare bit strings generated with the vec command rather than character strings. You could convert each character in the original string with a for loop and a hash table ie something like
      my $char_string = "AAATAT"; my $bit_string = ""; my %lookup_table = ("A" => "00", "T" => "01", "C" => "10", "G" => "11 +"); my @char_array=(split //, $char_string); for (my $i=0; $i <= $#char_array ; $i++) { vec($bit_string,$i,2) = $lookup_table{$char_array[$i]}; }

        There are several problem with that.

        1. Unless you were comparing each sequence against many others, the cost of encoding the sequences is going to exact a high price.
        2. Packing 4 chars per byte is a nice way of doing 1/4 of the work at the xor stage, if you can amortise the conversion. But DNA sequences often contain other characters, N and X for example, and I'm sure I've seen others. I'm not sure what they mean, but I've seen them.

          Once you are down to packing two, 3-bit or 4-bit values per byte, the cost of the encoding gets harder to claw back.

        3. And when you've done the xor, you have to count the results. With the string version, where the null byte means the characters were the same, and anything else, different, the "count the stars" mode of tr/// is an efficient way of achieving this.

          With the bit-encode sequences, you would have to count the zero and non-zero-ness of pairs, or triples or quads of bits. For 2 and 4 bits, this could be done using vec, but not very efficiently. Doing it for triples would be very laborious.


        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
        "Science is about questioning the status quo. Questioning authority".
        The "good enough" maybe good enough for the now, and perfection maybe unobtainable, but that should not preclude us from striving for perfection, when time, circumstance or desire allow.
      Or a little more plainly (update to explain: using the names the OP did, so he can see what corresponds):
      sub hd { my ($k, $l) = @_; my $diff = $k ^ $l; my $num_mismatch = $diff =~ tr/\0//c; }

      Caution: Contents may have been coded under pressure.

        I'm not sure that $k, $l are any more meaningful than $_[0], $_[1] and you pay a small performance penalty for obtaining those names.

        Does $diff really capture what it is naming? I've used $xor or $mask for that in the past, but I wonder if it isn't best left unnamed.

        I don't see any purpose in naming the return value, just to return it. Better to say return <EXPR>; and give the function a proper name like hammingDistance() and allow it to name the return value.

        While using tr/\0//c produces the same result as length - tr/\0/\0/, the major difference is that the former modifies the string being inspected, deleting the chars counted, where the latter does not. It just counts.

        On the short strings in the OP this is insignificant. But the technique works for any length strings, and as DNA strings can get very large indeed, the difference then becomes very significant.

        Overall, I prefer my version to yours, but each to their preference :)


        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
        "Science is about questioning the status quo. Questioning authority".
        The "good enough" maybe good enough for the now, and perfection maybe unobtainable, but that should not preclude us from striving for perfection, when time, circumstance or desire allow.
Re: Hamming Distance Between 2 Strings - Fast(est) Way?
by wazoox (Prior) on Oct 14, 2005 at 15:34 UTC

    Depending upon your data, simply checking first that the strings don't match can make it much faster :

    sub hd1 { #String length is assumed to be equal my ($k,$l) = @_; # next line added return 0 if $k eq $l; my $len = length ($k); my $num_mismatch = 0; for (my $i=0; $i<$len; $i++) { ++$num_mismatch if substr($k, $i, 1) ne substr($l, $i, 1); } return $num_mismatch; }
Re: Hamming Distance Between 2 Strings - Fast(est) Way?
by inman (Curate) on Oct 18, 2005 at 14:15 UTC
    This code uses an XOR to compare the two strings. This is case sensitive and works when the two strings are the same length.
    sub hd { return ($_[0] ^ $_[1]) =~ tr/\001-\255//; }
      Hello I'm confused with the benchmark by Edward. The XOR code seems to be 3 times faster than the 'mine' using the following benchmark (shouldn't benchmark code always be published?)
      #!/usr/bin/perl -w use strict; my $s1 = 'AAAAA'; my $s2 = 'ATCAA'; for (my $i=0;$i<600000;$i++){ #choose one of the two methods #hd($s1,$s2); # real 0m1.401s hd2($s1,$s2); # real 0m0.405s } sub hd{ my ($k,$l) = @_; my $len = length ($k); my $num_mismatch = 0; for (my $i=0; $i<$len; $i++) { ++$num_mismatch if substr($k, $i, 1) ne substr($l, $i, 1); } return $num_mismatch; } sub hd2 { return ($_[0] ^ $_[1]) =~ tr/\001-\255//; }
        Duh! I misread Edwards' benchmark. missed that the benchmark was /s (instead of s) - shouldn't I be reading more carefully? :-)
        Sorry this is probably a really dumb question...but.. what is the ^ character doing in this subfunction? I know in regex it means match the first... and its used as a way to negate character class when between brackets...but I cant find its use like shown above.

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: perlquestion [id://500235]
Approved by chester
Front-paged by kwaping
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others exploiting the Monastery: (5)
As of 2014-07-26 20:46 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    My favorite superfluous repetitious redundant duplicative phrase is:









    Results (179 votes), past polls