blast_kappa.c

ncbi源码

 
  Blast_HSPInit(queryStart, queryStart + queryAlignmentExtent,
                matchStart, matchStart + matchAlignmentExtent,
                0, 0, 0, 0, newScore, &editBlock, &(newSW->hsp));
  newSW->hsp->evalue = newEvalue;
 
  SWAlignGetNumIdentical(newSW, query_seq); /* Calculate num identities, attach to HSP. */

  self->alignments  = newSW;
}


/**
 * The struct SWheapRecord data type is used below to define the
 * internal structure of a SWheap (see below).  A SWheapRecord
 * represents all alignments of a query sequence to a particular
 * matching sequence.
 *
 * The SWResults::theseAlignments field is a linked list of alignments
 * of the query-subject pair.  The list is ordered by evalue in
 * descending order. Thus the first element has biggest (worst) evalue
 * and the last element has smallest (best) evalue. 
 */
typedef struct SWheapRecord {
  double bestEvalue;       /* best (smallest) evalue of all alignments
                                 * in the record */
  SWResults *theseAlignments;   /* a list of alignments */
} SWheapRecord;


/**  Compare two records in the heap.  
 * @param place1 the first record to be compared [in]
 * @param place2 the other record to be compared [in]
 */
static Boolean
SWheapRecordCompare(SWheapRecord * place1,
                    SWheapRecord * place2)
{
  return ((place1->bestEvalue > place2->bestEvalue) ||
          (place1->bestEvalue == place2->bestEvalue &&
           place1->theseAlignments->subject_index >
           place2->theseAlignments->subject_index));
}


/**  swap two records in the heap
 * @param heapArray holds the records to be swapped [in][out]
 * @param i the first record to be swapped [in]
 * @param j the other record to be swapped [in]
 */
static void
SWheapRecordSwap(SWheapRecord * heapArray,
                 Int4 i,
                 Int4 j)
{
  /* bestEvalue and theseAlignments are temporary variables used to
   * perform the swap. */
  double bestEvalue       = heapArray[i].bestEvalue;
  SWResults *theseAlignments   = heapArray[i].theseAlignments;

  heapArray[i].bestEvalue      = heapArray[j].bestEvalue;
  heapArray[i].theseAlignments = heapArray[j].theseAlignments;

  heapArray[j].bestEvalue      = bestEvalue;
  heapArray[j].theseAlignments = theseAlignments;
}


#ifdef KAPPA_INTENSE_DEBUG

/**
 * Verifies that the array heapArray[i] .. heapArray[n] is ordered so
 * as to be a valid heap.  This routine checks every element in the array,
 * an so is very time consuming.  It is for debugging purposes only.
 */
static Boolean
SWheapIsValid(SWheapRecord * heapArray,
              Int4 i,
              Int4 n)
{
  /* indices of nodes to the left and right of node i */
  Int4 left = 2 * i, right = 2 * i + 1;        

  if(right <= n) {
    return !SWheapRecordCompare(&(heapArray[right]), &(heapArray[i])) &&
      SWheapIsValid(heapArray, right, n);
  }
  if(left <= n) {
    return !SWheapRecordCompare(&(heapArray[left]), &(heapArray[i])) &&
      SWheapIsValid(heapArray, left, n);
  }
  return TRUE;
}

#define KAPPA_ASSERT(expr) ((expr) ? 0 : \
(fprintf( stderr, "KAPPA_ASSERT failed line %d: %s", __LINE__, #expr ), \
exit(1)))
#else
#define KAPPA_ASSERT(expr) (void)(0)
#endif


/** On entry, all but the first element of the array heapArray[i]
 * .. heapArray[n] are in valid heap order.  This routine rearranges
 * the elements so that on exit they all are in heap order. 
 * @param heapArray holds the heap [in][out]
 * @param i ?? [in]
 * @param n ?? [in]
 */
static void
SWheapifyDown(SWheapRecord * heapArray,
              Int4 i,
              Int4 n)
{
  Boolean moreswap = TRUE;      /* is more swapping needed */
  Int4 left, right, largest;    /* placeholders for indices in swapping */
  do {
    left  = 2 * i;
    right = 2 * i + 1;
    if((left <= n) &&
       (SWheapRecordCompare(&(heapArray[left]), &(heapArray[i]))))
      largest = left;
    else
      largest = i;
    if((right <= n) &&
       (SWheapRecordCompare(&(heapArray[right]), &(heapArray[largest]))))
      largest  = right;
    if(largest != i) {
      SWheapRecordSwap(heapArray, i, largest);
      /* push largest up the heap */
      i       = largest;       /* check next level down */
    } else
      moreswap = FALSE;
  } while(moreswap);            /* function builds the heap */
  KAPPA_ASSERT(SWheapIsValid(heapArray, i, n));
}


/** On entry, all but the last element of the array heapArray[i]
 *   .. heapArray[n] are in valid heap order.  This routine rearranges
 *   the elements so that on exit they all are in heap order.
 * @param heapArray holds the heap [in][out]
 * @param i the largest element to work with [in]
 * @param n the largest element in the heap [in]
 */
static void
SWheapifyUp(SWheapRecord * heapArray,
            Int4 i,
            Int4 n)
{
  Int4 parent = i / 2;          /* index to the node that is the
                                   parent of node i */
  while(parent >= 1 &&
        SWheapRecordCompare(&(heapArray[i]), &(heapArray[parent]))){
    SWheapRecordSwap(heapArray, i, parent);

    i       = parent;
    parent /= 2;
  }
  KAPPA_ASSERT(SWheapIsValid(heapArray, 1, n));
}

/** A SWheap represents a collection of alignments between one query
 * sequence and several matching subject sequences.  
 *
 * Each matching sequence is allocated one record in a SWheap.  The
 * eValue of a query-subject pair is the best (smallest positive)
 * evalue of all alignments between the two sequences.
 * 
 * A match will be inserted in the the SWheap if:
 * - there are fewer that SWheap::heapThreshold elements in the SWheap;
 * - the eValue of the match is <= SWheap::ecutoff; or
 * - the eValue of the match is less than the largest (worst) eValue
 *   already in the SWheap.
 *
 * If there are >= SWheap::heapThreshold matches already in the SWheap
 * when a new match is to be inserted, then the match with the largest
 * (worst) eValue is removed, unless the largest eValue <=
 * SWheap::ecutoff.  Matches with eValue <= SWheap::ecutoff are never
 * removed by the insertion routine.  As a consequence, the SWheap can
 * hold an arbitrarily large number of matches, although it is
 * atypical for the number of matches to be greater than
 * SWheap::heapThreshold.
 *
 * Once all matches have been collected, the SWheapToFlatList routine
 * may be invoked to return a list of all alignments. (see below).
 *
 * While the number of elements in a heap < SWheap::heapThreshold, the
 * SWheap is implemented as an unordered array, rather than a
 * heap-ordered array.  The SWheap is converted to a heap-ordered
 * array as soon as it becomes necessary to order the matches by
 * evalue.  The routines that operate on a SWheap should behave
 * properly whichever state the SWheap is in.
 */
struct SWheap {
  Int4 n;                       /**< The current number of elements */
  Int4 capacity;                /**< The maximum number of elements that may be 
                                   inserted before the SWheap must be resized */
  Int4 heapThreshold;           /**< see above */
  double ecutoff;          /**< matches with evalue below ecutoff may
                                   always be inserted in the SWheap */
  double worstEvalue;      /**< the worst (biggest) evalue currently in
                                   the heap */

  SWheapRecord *array;          /**< the SWheapRecord array if the SWheap is
                                   being represented as an unordered array */
  SWheapRecord *heapArray;      /**< the SWheapRecord array if the SWheap is
                                   being represented as an heap-ordered
                                   array. At least one of (array, heapArray)
                                   is NULL */

};
typedef struct SWheap SWheap;


/** Convert a SWheap from a representation as an unordered array to
 *    a representation as a heap-ordered array. 
 * @param self record to be modified [in][out]
 */
static void
ConvertToHeap(SWheap * self)
{
  if(NULL != self->array) {
    Int4 i;                     /* heap node index */
    Int4 n;                     /* number of elements in the heap */
    /* We aren't already a heap */
    self->heapArray = self->array;
    self->array     = NULL;

    n = self->n;
    for(i = n / 2; i >= 1; --i) {
      SWheapifyDown(self->heapArray, i, n);
    }
  }
  KAPPA_ASSERT(SWheapIsValid(self->heapArray, 1, self->n));
}

/*When the heap is about to exceed its capacity, it will be grown by
 *the minimum of a multiplicative factor of SWHEAP_RESIZE_FACTOR
 *and an additive factor of SWHEAP_MIN_RESIZE. The heap never
 *decreases in size */
#define SWHEAP_RESIZE_FACTOR 1.5
#define SWHEAP_MIN_RESIZE 100

/** Return true if self would insert a match that had the given eValue 
 * @param self record to be modified [in][out]
 * @param eValue specified expect value [in]
 */
static Boolean
SWheapWouldInsert(SWheap * self,
                  double eValue)
{
  return self->n < self->heapThreshold ||
    eValue <= self->ecutoff ||
    eValue < self->worstEvalue;
}


/** Try to insert matchRecord into the SWheap. The alignments stored in
 * matchRecord are used directly, i.e. they are not copied, but are
 * rather stored in the SWheap or deleted
 * @param self record to be modified [in][out]
 * @param matchRecord record to be inserted [in]
 */
static void
SWheapInsert(SWheap * self,
             Kappa_MatchRecord * matchRecord)
{
  if(self->array && self->n >= self->heapThreshold) {
    ConvertToHeap(self);
  }
  if(self->array != NULL) {
    /* "self" is currently a list. Add the new alignments to the end */
    SWheapRecord *heapRecord;   /* destination for the new alignments */ 
    heapRecord                  = &self->array[++self->n];
    heapRecord->bestEvalue      = matchRecord->eValue;
    heapRecord->theseAlignments = matchRecord->alignments;
    if( self->worstEvalue < matchRecord->eValue ) {
      self->worstEvalue = matchRecord->eValue;
    }
  } else {                      /* "self" is currently a heap */
    if(self->n < self->heapThreshold ||
       (matchRecord->eValue <= self->ecutoff &&
        self->worstEvalue <= self->ecutoff)) {
      SWheapRecord *heapRecord; /* Destination for the new alignments */
      /* The new alignments must be inserted into the heap, and all old
       * alignments retained */
      if(self->n >= self->capacity) {
        /* The heap must be resized */
        Int4 newCapacity;       /* capacity the heap will have after
                                 * it is resized */
        newCapacity      = MAX(SWHEAP_MIN_RESIZE + self->capacity,
                               (Int4) (SWHEAP_RESIZE_FACTOR * self->capacity));
        self->heapArray  = (SWheapRecord *)
          realloc(self->heapArray, (newCapacity + 1) * sizeof(SWheapRecord));
        self->capacity   = newCapacity;
      }
      /* end if the heap must be resized */
      heapRecord    = &self->heapArray[++self->n];
      heapRecord->bestEvalue      = matchRecord->eValue;
      heapRecord->theseAlignments = matchRecord->alignments;

      SWheapifyUp(self->heapArray, self->n, self->n);
    } else {
      /* Some set of alignments must be discarded */
      SWResults *discardedAlignments = NULL;      /* alignments that
                                                   * will be discarded
                                                   * so that the new
                                                   * alignments may be
                                                   * inserted. */

      if(matchRecord->eValue >= self->worstEvalue) {
        /* the new alignments must be discarded */
        discardedAlignments = matchRecord->alignments;
      } else {
        /* the largest element in the heap must be discarded */
        SWheapRecord *heapRecord;     /* destination for the new alignments */
        discardedAlignments         = self->heapArray[1].theseAlignments;

        heapRecord                  = &self->heapArray[1];
        heapRecord->bestEvalue      = matchRecord->eValue;
        heapRecord->theseAlignments = matchRecord->alignments;

        SWheapifyDown(self->heapArray, 1, self->n);
      }
      /* end else the largest element in the heap must be discarded */
      while(discardedAlignments != NULL) {
        /* There are discarded alignments that have not been freed */
        SWResults *thisAlignment;     /* the head of the list of
                                       * discarded alignments */
        thisAlignment        = discardedAlignments;
        discardedAlignments  = thisAlignment->next;
        sfree(thisAlignment);
      }
      /* end while there are discarded alignments that have not been freed */
    } 
    /* end else some set of alignments must be discarded */