lookup_util.c
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:5k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: lookup_util.c,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 18:08:08  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.8
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /* $Id: lookup_util.c,v 1000.1 2004/06/01 18:08:08 gouriano Exp $
  10.  * ===========================================================================
  11.  *
  12.  *                            PUBLIC DOMAIN NOTICE
  13.  *               National Center for Biotechnology Information
  14.  *
  15.  *  This software/database is a "United States Government Work" under the
  16.  *  terms of the United States Copyright Act.  It was written as part of
  17.  *  the author's offical duties as a United States Government employee and
  18.  *  thus cannot be copyrighted.  This software/database is freely available
  19.  *  to the public for use. The National Library of Medicine and the U.S.
  20.  *  Government have not placed any restriction on its use or reproduction.
  21.  *
  22.  *  Although all reasonable efforts have been taken to ensure the accuracy
  23.  *  and reliability of the software and data, the NLM and the U.S.
  24.  *  Government do not and cannot warrant the performance or results that
  25.  *  may be obtained by using this software or data. The NLM and the U.S.
  26.  *  Government disclaim all warranties, express or implied, including
  27.  *  warranties of performance, merchantability or fitness for any particular
  28.  *  purpose.
  29.  *
  30.  *  Please cite the author in any work or product based on this material.
  31.  *
  32.  * ===========================================================================
  33.  */
  34. /** @file lookup_util.c
  35.  * @todo FIXME needs file description
  36.  */
  37. static char const rcsid[] = 
  38.     "$Id: lookup_util.c,v 1000.1 2004/06/01 18:08:08 gouriano Exp $";
  39. #include <algo/blast/core/lookup_util.h>
  40. void __sfree(void **x)
  41. {
  42.   free(*x);
  43.   *x=NULL;
  44.   return;
  45. }
  46. static void fkm_output(Int4* a, Int4 n, Int4 p, Uint1* output, Int4* cursor, Uint1* alphabet);
  47. static void fkm(Int4* a, Int4 n, Int4 k, Uint1* output, Int4* cursor, Uint1* alphabet);
  48. Int4 iexp(Int4 x, Int4 n)
  49. {
  50. Int4 r,y;
  51. r=1;
  52. y=x;
  53.  if(n==0) return 1;
  54.  if(x==0) return 0;
  55.   
  56. while (n > 1)
  57.   { 
  58.     if ( (n%2)==1) r *= y;
  59.     n = n >> 1;
  60.     y = y*y;
  61.   }
  62.  r = r*y;
  63.  return r;
  64. }
  65. Int4 ilog2(Int4 x)
  66. {
  67.   Int4 lg=0;
  68.   if (x==0) return 0;
  69.   while ( ( x = x >> 1) )
  70.     lg++;
  71.   return lg;
  72. }
  73. Int4 makemask(Int4 x)
  74. {
  75. Int4 mask=1;
  76. if (x==0) return 0;
  77. mask = mask << (x-1);
  78. mask |= (mask >> 1);
  79. mask |= (mask >> 2);
  80. mask |= (mask >> 4);
  81. mask |= (mask >> 8);
  82. mask |= (mask >> 16);
  83. return mask;
  84. }
  85. /** Output a Lyndon word as part of a de Bruijn sequence.
  86.  *
  87.  *if the length of a lyndon word is divisible by n, print it. 
  88.  * @param a the shift register
  89.  * @param p
  90.  * @param n
  91.  * @param output the output sequence
  92.  * @param cursor current location in the output sequence
  93.  * @param alphabet optional translation alphabet
  94.  */
  95. static void fkm_output(Int4* a, Int4 n, Int4 p, Uint1* output, Int4* cursor, Uint1* alphabet)
  96. {
  97.   Int4 i;
  98.   if (n % p == 0)
  99.     for(i=1;i<=p;i++)
  100.       {
  101. if (alphabet != NULL)
  102.   output[*cursor] = alphabet[ a[i] ];
  103. else
  104.   output[*cursor] = a[i];
  105. *cursor = *cursor + 1;
  106.       }
  107. }
  108. /**
  109.  * iterative fredricksen-kessler-maiorana algorithm
  110.  * to generate de bruijn sequences.
  111.  *
  112.  * generates all lyndon words, in lexicographic order.
  113.  * the concatenation of all lyndon words whose length is 
  114.  * divisible by n yields a de bruijn sequence.
  115.  *
  116.  * further, the sequence generated is of shortest lexicographic length.
  117.  *
  118.  * references: 
  119.  * http://mathworld.wolfram.com/deBruijnSequence.html
  120.  * http://www.theory.csc.uvic.ca/~cos/inf/neck/NecklaceInfo.html
  121.  * http://www.cs.usyd.edu.au/~algo4301/ , chapter 7.2
  122.  * http://citeseer.nj.nec.com/ruskey92generating.html
  123.  *
  124.  * @param a the shift register
  125.  * @param n the number of letters in each word
  126.  * @param k the size of the alphabet
  127.  * @param output the output sequence
  128.  * @param cursor the current location in the output sequence
  129.  * @param alphabet optional translation alphabet
  130.  */
  131. static void fkm(Int4* a, Int4 n, Int4 k, Uint1* output, Int4* cursor, Uint1* alphabet)
  132. {
  133.   Int4 i,j;
  134.   fkm_output(a,n,1,output,cursor,alphabet);
  135.   i=n;
  136. do
  137.   {
  138.     a[i] = a[i] + 1;
  139.     for(j=1;j<=n-i;j++)
  140.       a[j+i] = a[j];
  141.     fkm_output(a,n,i,output,cursor,alphabet);
  142.     i=n;
  143.     while (a[i] == k-1)
  144.       i--;
  145.   }
  146.  while (i>0);
  147. }
  148. void debruijn(Int4 n, Int4 k, Uint1* output, Uint1* alphabet)
  149. {
  150.   Int4* a;
  151.   Int4 cursor=0;
  152.   /* n+1 because the array is indexed from one, not zero */
  153.   a = (Int4*) calloc( (n+1), sizeof(Int4));
  154.   /* compute the (n,k) de Bruijn sequence and store it in output */
  155.  
  156.   fkm(a,n,k,output,&cursor,alphabet);
  157.   
  158.   sfree(a);
  159.   return;
  160. }