website/archive/binaries/mac/src/clustalo/src/squid/rk.c

   1 /*****************************************************************
   2  * SQUID - a library of functions for biological sequence analysis
   3  * Copyright (C) 1992-2002 Washington University School of Medicine
   4  *
   5  *     This source code is freely distributed under the terms of the
   6  *     GNU General Public License. See the files COPYRIGHT and LICENSE
   7  *     for details.
   8  *****************************************************************/
   9
  10 /* rk.c (originally from rnabob's patsearch.c)
  11  *
  12  * Contains a compiler and a search engine for Rabin-Karp
  13  * based primary sequence pattern searching on encoded
  14  * sequences.
  15  *
  16  * See Sedgewick, _Algorithms_, for a general discussion of
  17  * the Rabin-Karp algorithm. See the rkcomp or rkexec man
  18  * pages for specific details.
  19  *
  20  * RCS $Id: rk.c 217 2011-03-19 10:27:10Z andreas $ (Original squid RCS Id: rk.c,v 1.2 1998/10/09 18:07:16 eddy Exp)
  21  */
  22
  23 #include <stdio.h>
  24 #include <string.h>
  25 #include <ctype.h>
  26 #include "squid.h"              /* seq encoding utilities and typedefs */
  27 #include "rk.h"
  28
  29
  30 #ifdef MEMDEBUG
  31 #include "dbmalloc.h"
  32 #endif
  33
  34 Hashseq
  35 rkcomp(char *probe)               /* A,C,G,T/U, N probe string, 0-8 nt long */
  36 {
  37   Hashseq hashprobe = 0;
  38   char    coded[RK_HASHSIZE + 1];
  39   int     len;
  40   int     i;
  41                                 /* check bounds violation on probe */
  42   if ((len =  strlen(probe)) > RK_HASHSIZE) return 0;
  43                                 /* encode the probe */
  44   if (seqencode(coded, probe) == 0) return 0;
  45                                 /* pack the probe into a Hashseq */
  46   for (i = 0; i < len; i++)
  47     {
  48       hashprobe <<= 4;
  49       hashprobe |= (Hashseq) coded[i];
  50     }
  51                                 /* left adjust as needed */
  52   for (; i < RK_HASHSIZE; i++)
  53     {
  54       hashprobe <<= 4;
  55       hashprobe |= (Hashseq) NTN;
  56     }
  57                                 /* return the compiled probe */
  58   return hashprobe;
  59 }
  60
  61 int
  62 rkseq(Hashseq   hashprobe,      /* up to 8 nt packed into the probe */
  63       char     *sequence)       /* encoded sequence                 */
  64 {
  65   long     i;
  66   long     pos = 0;
  67   Hashseq  target = 0;
  68
  69                                 /* initialize the target hashseq */
  70   for (i = 0; i < RK_HASHSIZE; i++)
  71     {
  72       if (*(sequence + i) == NTEND)
  73         break;
  74       target <<= 4;
  75       target |=  (Hashseq) (*(sequence + i));
  76     }
  77
  78   while (*(sequence + pos + RK_HASHSIZE -1) != NTEND)
  79     {
  80 #ifdef DEBUG
  81       printf("hashprobe: ");
  82       writehash(hashprobe);
  83       printf("\ttarget: ");
  84       writehash(target);
  85       printf("\nhashprobe & target: ");
  86       writehash(hashprobe & target);
  87       printf("\n");
  88 #endif
  89       if ((hashprobe & target) == target)
  90         return ((int) pos);
  91       target <<= 4;
  92       target |=  (Hashseq) (*(sequence + pos + RK_HASHSIZE));
  93       pos++;
  94     }
  95                                 /* now we deal with an end effect */
  96   for (i = 0; i < RK_HASHSIZE; i++)
  97     {
  98       target |= (Hashseq) NTN;
  99       if ((hashprobe & target) == target)
 100         return ((int) pos);
 101       target <<=4;
 102       pos++;
 103     }
 104
 105   return(-1);
 106 }
 107
 108
 109 #ifdef DEBUG                    /* Debugging aids */
 110
 111 static void
 112 writehash(Hashseq   hashseq)
 113 {
 114   int idx;
 115   int sym;
 116
 117   if (hashseq/16)
 118     writehash(hashseq/16);
 119
 120   sym = (int) (hashseq % 16);
 121   if (sym == 0)
 122     putchar('-');
 123   else
 124     {
 125       for (idx = 0; sym != iupac[idx].code && idx < IUPACSYMNUM; idx++);
 126       if (idx > IUPACSYMNUM)
 127         printf("(%d)", sym);
 128       else
 129         putchar(iupac[idx].sym);
 130     }
 131 }
 132
 133 #endif