code_perl/Maasha/Backtrack.pm

   1 package Maasha::Backtrack;
   2
   3 # Copyright (C) 2006-2011 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # Routines for locating patterns in sequences.
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use warnings;
  32 use strict;
  33 use Exporter;
  34 use Data::Dumper;
  35 use Maasha::Common;
  36 use Maasha::Seq;
  37 use Maasha::Fasta;
  38 use vars qw ( @ISA @EXPORT );
  39
  40 @ISA = qw( Exporter );
  41
  42 use constant {
  43     SEQ_NAME => 0,
  44     SEQ  => 1,
  45 };
  46
  47
  48 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  49
  50
  51 use Inline (C => <<'END_C');
  52 #define MATCH(A,B) ((bitmap[(unsigned short) A] & bitmap[(unsigned short) B]) != 0)
  53
  54 // IUPAC ambiguity codes
  55 unsigned short bitmap[256] = {
  56     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  57     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  58     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  59     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  60     0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
  61     0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
  62     0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
  63     0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
  64     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  65     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  66     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  67     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  68     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  69     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  70     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  71     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
  72 };
  73
  74 // ss is the start of the string, used only for reporting the match endpoints.
  75 unsigned int backtrack(char* ss, char* s, char* p, int mm, int ins, int del)
  76 {
  77     unsigned int r = 0;
  78
  79     while (*s && MATCH(*s, *p)) ++s, ++p;    // OK to always match longest segment
  80
  81     if (!*p)
  82         return (unsigned int) (s - ss) - 1;
  83     else
  84     {
  85         if (mm && *s && *p && (r = backtrack(ss, s + 1, p + 1, mm - 1, ins, del))) return r;
  86         if (ins && *s &&      (r = backtrack(ss, s + 1, p, mm, ins - 1, del)))     return r;
  87         if (del && *p &&      (r = backtrack(ss, s, p + 1, mm, ins, del - 1)))     return r;
  88     }
  89
  90     return 0;
  91 }
  92
  93 // Find match of p in s, with at most mm mismatches, ins insertions and del deletions.
  94 void patscan(char* s, char* p, int pos, int mm, int ins, int del)
  95 {
  96     char*        ss = s;
  97     unsigned int e;
  98
  99     while (*s && pos) ++s, --pos;   // Fast forward until pos
 100
 101     while (*s)
 102     {
 103         e = backtrack(ss, s, p, mm, ins, del);
 104
 105         if (e)
 106         {
 107             // printf("beg: %d   end: %d\n", (int) (s - ss), e);
 108             Inline_Stack_Vars;
 109             Inline_Stack_Reset;
 110             Inline_Stack_Push( sv_2mortal( newSViv(s - ss)));
 111             Inline_Stack_Push( sv_2mortal( newSViv(e)));
 112             Inline_Stack_Done;
 113         }
 114
 115         s++;
 116     }
 117 }
 118
 119 END_C
 120
 121
 122 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<