Beefy Boxes and Bandwidth Generously Provided by pair Networks Bob
Problems? Is your data what you think it is?
 
PerlMonks  

Re^2: searching for a pattern using suffix arrays (close)

by abualiga (Scribe)
on Oct 09, 2013 at 02:19 UTC ( #1057473=note: print w/ replies, xml ) Need Help??


in reply to Re: searching for a pattern using suffix arrays (close)
in thread searching for a pattern using suffix arrays

Greatly appreciate your input Tye!

I accepted your corrections, and will adopt some of your style practices. I messed with the binary search sub a bit and may have stumbled upon a way to get the first index. I tried several random DNA strings, and it works. Please see below and let me know your thoughts, at your convenience.

thanks!

#!/usr/bin/perl use Modern::Perl '2011'; use autodie; my $str = 'atgagatacgattacagattacagatcatacgataccatacgatc$'; my @suff; # hold suffixes while ( 1 < length $str ) { # condition excludes '$' from @suff push @suff, $str; substr ( $str, 0, 1, '' ); } # lexically ordered suffix array my @indices = sort { $suff[$a] cmp $suff[$b] } 0 .. $#suff; for ( 0 .. $#indices ) { say $indices[$_]+1, ": ", $suff[ $indices[$_] ]; } my $pattern = 'tac'; # find smallest index where pattern matches first n characters of suff +ix my $start = bsearch( \@indices, $pattern ); # get consecutive positions where pattern matches first n chars of suf +fix my @positions; for my $index ( $start .. $#indices ) { if ( $suff[ $indices[ $index ] ] =~ /^$pattern/ ) { push @positions, $indices[ $index ] + 1; } else { last; } } # print results say "\nPattern \"$pattern\" found at (position: suffix):"; for my $pos ( sort{ $a <=> $b } @positions ) { chop( my $s = $suff[ $pos - 1] ); say "$pos: $s"; } sub bsearch { my ( $indref, $pat ) = @_; my $mid; my ( $lo, $hi ) = ( 0, $#$indref ); while ( 1 ) { $mid = int( ( $lo + $hi ) / 2 ); if ( $hi == $lo ) { return $mid } if( ( $pattern cmp $suff[ $indices[ $mid ] ] ) < 0 ) { $hi = $mid; } else { $lo = $mid + 1; } } return "pattern not found in string!"; } #output abualiga:~$ ./suffixArray.pl 21: acagatcatacgataccatacgatc$ 14: acagattacagatcatacgataccatacgatc$ 35: accatacgatc$ 30: acgataccatacgatc$ 40: acgatc$ 8: acgattacagattacagatcatacgataccatacgatc$ 4: agatacgattacagattacagatcatacgataccatacgatc$ 23: agatcatacgataccatacgatc$ 16: agattacagatcatacgataccatacgatc$ 33: ataccatacgatc$ 28: atacgataccatacgatc$ 38: atacgatc$ 6: atacgattacagattacagatcatacgataccatacgatc$ 43: atc$ 25: atcatacgataccatacgatc$ 1: atgagatacgattacagattacagatcatacgataccatacgatc$ 18: attacagatcatacgataccatacgatc$ 11: attacagattacagatcatacgataccatacgatc$ 45: c$ 22: cagatcatacgataccatacgatc$ 15: cagattacagatcatacgataccatacgatc$ 27: catacgataccatacgatc$ 37: catacgatc$ 36: ccatacgatc$ 31: cgataccatacgatc$ 41: cgatc$ 9: cgattacagattacagatcatacgataccatacgatc$ 3: gagatacgattacagattacagatcatacgataccatacgatc$ 32: gataccatacgatc$ 5: gatacgattacagattacagatcatacgataccatacgatc$ 42: gatc$ 24: gatcatacgataccatacgatc$ 17: gattacagatcatacgataccatacgatc$ 10: gattacagattacagatcatacgataccatacgatc$ 20: tacagatcatacgataccatacgatc$ 13: tacagattacagatcatacgataccatacgatc$ 34: taccatacgatc$ 29: tacgataccatacgatc$ 39: tacgatc$ 7: tacgattacagattacagatcatacgataccatacgatc$ 44: tc$ 26: tcatacgataccatacgatc$ 2: tgagatacgattacagattacagatcatacgataccatacgatc$ 19: ttacagatcatacgataccatacgatc$ 12: ttacagattacagatcatacgataccatacgatc$ Pattern "tac" found at (position: suffix): 7: tacgattacagattacagatcatacgataccatacgatc 13: tacagattacagatcatacgataccatacgatc 20: tacagatcatacgataccatacgatc 29: tacgataccatacgatc 34: taccatacgatc 39: tacgatc


Comment on Re^2: searching for a pattern using suffix arrays (close)
Download Code
Re^3: searching for a pattern using suffix arrays (better)
by tye (Cardinal) on Oct 09, 2013 at 02:46 UTC

    Yeah, I can see value in getting the binary search to find the "first" potential match; especially how you've avoided using the regex within it.

    Many decent improvements.

    bsearch() is still accessing @suff like a global variable. You may have decided that that is fine.

    I'd rewrite:

    $mid = int( ( $lo + $hi ) / 2 ); if ( $hi == $lo ) { return $mid }

    as

    return $hi if $hi == $lo; $mid = int( ( $lo + $hi ) / 2 );

    though not for any compelling reasons.

    Your bsearch() now returns the first potential match (the first suffix after or equal to your search), so you might want to adjust your output for the case of no matches.

    The return value from bsearch() might be an error message but you unconditionally use it as a number. If Modern::Perl doesn't turn on warnings, then you should probably do that.

    if ( $suff[ $indices[ $index ] ] =~ /^$pattern/ ) { push @positions, $indices[ $index ] + 1; } else { last; }

    can be shortened as:

    last if $suff[ $indices[ $index ] ] !~ /^$pattern/; push @positions, $indices[ $index ] + 1;

    - tye        

      Thanks again Tye!

      I incorporated your changes: bsearch() is no longer accessing global variables, output is adjusted for a 'no match' case, and Modern::Perl does turn on warnings by default.

      #!/usr/bin/perl use Modern::Perl '2011'; use autodie; # pattern searching using suffix arrays and binary search open my $fh, '<', shift; chomp( my( $pattern, $str ) = <$fh> ); $str .= '$'; my @suff; # hold suffixes while ( 1 < length $str ) { # condition excludes '$' from @suff push @suff, $str; substr ( $str, 0, 1, '' ); } # lexically ordered suffix array my @indices = sort { $suff[$a] cmp $suff[$b] } 0 .. $#suff; #for ( 0 .. $#indices ) { # say $indices[$_]+1, ": ", $suff[ $indices[$_] ]; #} my $start = bsearch( \ @indices, $pattern, \ @suff ); # get consecutive positions, if any, # where pattern matches first n chars of suffix. my @positions; for my $index ( $start .. $#indices ) { last if $suff[ $indices[ $index ] ] !~ /^$pattern/; push @positions, $indices[ $index ] + 1; # omit +1 if 0-based inde +xing } # print results if ( @positions ) { # check if a match exists at all say "\nPattern \"$pattern\" found at position:"; for my $pos ( sort{ $a <=> $b } @positions ) { #chop( my $s = $suff[ $pos - 1] ); #suffix print "$pos "; } } else { say "\nPattern \"$pattern\" not found in string"; } print "\n"; # binary search. # find first potential match, ie first suffix after or equal to patter +n, # such that pattern potentially matches first n characters of suffix. sub bsearch { my ( $indref, $pat, $sufref ) = @_; my $mid; my ( $lo, $hi ) = ( 0, $#$indref ); while ( 1 ) { $mid = int( ( $lo + $hi ) / 2 ); return $mid if $hi == $lo; if( ( $pattern cmp $$sufref[ $$indref[ $mid ] ] ) < 0 ) { $hi = $mid; } else { $lo = $mid + 1; } } }

      I did a simple benchmark relative to what you would call a 'naive' pattern search, after checking that both give the same results, with a 9-char pattern and 8569-char string. I like the speed up.

      #!/usr/bin/perl use Modern::Perl '2011'; use autodie; use Benchmark qw/ timethese /; # benchmark pattern searching using naive approach vs suffix arrays open my $fh, '<', shift; chomp( my( $pattern, $str ) = <$fh> ); my $naive = sub { my @positions; while ( $str =~ /(?=($pattern))/g ) { push @positions, $-[0];; } #say "@positions"; }; my $sufbin = sub { $str .= '$'; my @suff; # hold suffixes while ( 1 < length $str ) { # condition excludes '$' from @suff push @suff, $str; substr ( $str, 0, 1, '' ); } # lexically ordered suffix array my @indices = sort { $suff[$a] cmp $suff[$b] } 0 .. $#suff; #for ( 0 .. $#indices ) { # say $indices[$_]+1, ": ", $suff[ $indices[$_] ]; #} my $start = bsearch( \ @indices, $pattern, \ @suff ); # get consecutive positions, if any, # where pattern matches first n chars of suffix. my @positions; for my $index ( $start .. $#indices ) { last if $suff[ $indices[ $index ] ] !~ /^$pattern/; push @positions, $indices[ $index ];# + 1; # omit +1 if 0-base +d indexing } # binary search. # find first potential match, ie first suffix after or equal to pa +ttern, # such that pattern potentially matches first n characters of suff +ix. sub bsearch { my ( $indref, $pat, $sufref ) = @_; my $mid; my ( $lo, $hi ) = ( 0, $#$indref ); while ( 1 ) { $mid = int( ( $lo + $hi ) / 2 ); return $mid if $hi == $lo; if( ( $pattern cmp $$sufref[ $$indref[ $mid ] ] ) < 0 ) { $hi = $mid; } else { $lo = $mid + 1; } } } }; timethese( -5, { Suffixbinary => $sufbin, Naive => $naive, } ); #output abualiga:~$ ./benchmarkPatternSearch.pl patternSearchData.txt Benchmark: running Naive, Suffixbinary for at least 5 CPU seconds... Naive: 6 wallclock secs ( 5.32 usr + 0.00 sys = 5.32 CPU) @ 92 +5.75/s (n=4925) Suffixbinary: 6 wallclock secs ( 5.33 usr + 0.00 sys = 5.33 CPU) @ +232758.91/s (n=1240605)

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://1057473]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others exploiting the Monastery: (4)
As of 2014-04-20 13:38 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    April first is:







    Results (485 votes), past polls