added backtrack.rb

git-svn-id: http://biopieces.googlecode.com/svn/trunk@1640 74ccb610-7750-0410-82ae-013aeee3265d

diff --git a/code_ruby/lib/maasha/backtrack.rb b/code_ruby/lib/maasha/backtrack.rb
new file mode 100644 (file)
index 0000000..a124895
--- /dev/null
+++ b/code_ruby/lib/maasha/backtrack.rb
@@ -0,0 +1,156 @@
+# Copyright (C) 2007-2011 Martin A. Hansen.
+
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 2
+# of the License, or (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+# http://www.gnu.org/copyleft/gpl.html
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+# This software is part of the Biopieces framework (www.biopieces.org).
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+require 'inline'
+
+# Module containing code to locate nucleotide patterns in sequences allowing for
+# ambiguity codes and a given maximum mismatches, insertions, and deletions. The
+# pattern match engine is based on a backtrack algorithm.
+module BackTrack
+  # ------------------------------------------------------------------------------
+  #   str.scan(pattern[, max_mismatches [, max_insertions [,max_deletions]]])
+  #   -> Array
+  #   str.scan(pattern[, max_mismatches [, max_insertions [,max_deletions]]]) { |match|
+  #     block
+  #   }
+  #   -> Match
+  #
+  # ------------------------------------------------------------------------------
+  # Method to iterate through a sequence to locate pattern matches starting at a
+  # given offset allowing for a maximum number of mismatches, insertions, and
+  # deletions. Matches found in block context return the Match object. Otherwise
+  # matches are returned in an Array of Match objects.
+  def patscan(pattern, offset = 0, max_mismatches = 0, max_insertions = 0, max_deletions = 0)
+    matches = []
+
+    if block_given?
+      while result = self.track(pattern, offset, max_mismatches, max_insertions, max_deletions)
+        yield Match.new(result.first, result.last, self.seq[result.first ... result.first + result.last])
+        offset = result.first + 1
+      end
+    else
+      while result = self.track(pattern, offset, max_mismatches, max_insertions, max_deletions)
+        matches << Match.new(result.first, result.last, self.seq[result.first ... result.first + result.last])
+        offset = result.first + 1
+      end
+    end
+
+    matches unless block_given?
+  end
+
+  private
+
+  inline do |builder|
+    builder.add_static "id_seq", 'rb_intern("@seq")', "ID"
+
+    builder.prefix %{
+      #define MATCH(A,B) ((bitmap[(unsigned short) A] & bitmap[(unsigned short) B]) != 0)
+    }
+
+    builder.prefix %{
+      unsigned short bitmap[256] = {
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
+          0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
+          0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
+          0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+      };
+    }
+
+    # ss is the start of the string, used only for reporting the match endpoints.
+    builder.prefix %{
+      unsigned int backtrack(char* ss, char* s, char* p, int mm, int ins, int del)
+      {
+          unsigned int r = 0;
+
+          while (*s && MATCH(*s, *p)) ++s, ++p;    // OK to always match longest segment
+
+          if (!*p)
+              return (unsigned int) (s - ss);
+          else
+          {
+              if (mm && *s && *p && (r = backtrack(ss, s + 1, p + 1, mm - 1, ins, del))) return r;
+              if (ins && *s &&      (r = backtrack(ss, s + 1, p, mm, ins - 1, del)))     return r;
+              if (del && *p &&      (r = backtrack(ss, s, p + 1, mm, ins, del - 1)))     return r;
+          }
+
+          return 0;
+      }
+    }
+ 
+    # Find all occurrences of p starting at pos in s, with at most
+    # mm mismatches, ins insertions and del deletions.
+    builder.c %{
+      static VALUE track(char* p, int pos, int mm, int ins, int del)
+      {
+          VALUE        seq = rb_ivar_get(self, id_seq);
+          char*        s   = StringValuePtr(seq);
+          char*        ss  = s;
+          unsigned int e;
+          VALUE        tuple;
+
+          while (*s && pos) ++s, --pos;   // Fast forward until pos
+
+          while (*s)
+          {
+              e = backtrack(ss, s, p, mm, ins, del);
+
+              if (e)
+              {
+                  tuple = rb_ary_new();
+                  rb_ary_push(tuple, INT2FIX((int) (s - ss)));
+                  rb_ary_push(tuple, INT2FIX((int) e - (s - ss)));
+                  return tuple;
+              }
+
+              s++;
+          }
+       }
+    }
+  end
+
+  class Match
+    attr_reader :pos, :length, :match
+
+    def initialize(pos, length, match)
+      @pos    = pos
+      @length = length
+      @match  = match
+    end
+  end
+end
+
+
+__END__