a44421a79fb36cc2036fe116b97ea3bc9590cd0c
braney
  Fri Dec 2 09:34:39 2011 -0800
removed rcsid (#295)
diff --git src/lib/pairHmm.c src/lib/pairHmm.c
index e7c9464..67eda43 100644
--- src/lib/pairHmm.c
+++ src/lib/pairHmm.c
@@ -1,347 +1,346 @@
 /* pairHmm - stuff to help implement pairwise hidden markov models,
  * which are useful ways of aligning two sequences. 
  *
  * This file is copyright 2000-2004 Jim Kent, but license is hereby
  * granted for all use - public, private or commercial. */
 
 #include "common.h"
 #include "axt.h"
 #include "pairHmm.h"
 
-static char const rcsid[] = "$Id: pairHmm.c,v 1.5 2005/04/11 06:43:58 kent Exp $";
 
 UBYTE phmmNullMommy = 0; /* mommy value for orphans.... */
 
 void phmmUnpackMommy(UBYTE mommy, int *retStateIx, int *retQoff, 
 	int *retToff)
 /* Unpack state, query, and target offset. */
 {
 *retStateIx = (mommy&31);
 *retQoff = -((mommy&32)>>5);
 *retToff = -((mommy&64)>>6);
 }
 
 struct phmmMatrix *phmmMatrixNew(int stateCount,
     char *query, int querySize, char *target, int targetSize)
 /* Allocate all memory required for an phmmMatrix. Set up dimensions. */
 {
 int i;
 struct phmmMommy *allCells;
 int allCellSize;
 int rowSize;
 int *allScores;
 struct phmmMatrix *am;
 
 AllocVar(am);
 am->query = query;
 am->target = target;
 am->querySize = querySize;
 am->targetSize = targetSize;
 am->qDim = rowSize = am->querySize + 1;
 am->tDim = am->targetSize + 1;
 am->stateCount = stateCount;
 am->stateSize = rowSize * am->tDim;
 am->stateByteSize = am->stateSize * sizeof(struct phmmMommy);
 am->states = needMem(stateCount * sizeof(struct phmmState));
 
 /* Initialize matrix of cells for each state. */
 allCellSize = stateCount * am->stateByteSize;
 am->allCells = allCells = needLargeMem(allCellSize); 
 memset(allCells, 0, allCellSize);
 for (i=0; i<stateCount; ++i)
     {
     am->states[i].cells = allCells;
     allCells += am->stateSize;
     am->states[i].stateIx = i;
     }
 
 /* Initialize two rows of scores for each state. */
 allScores = am->allScores = needMem(rowSize * 2 * stateCount * sizeof(int) );
 for (i=0; i<stateCount; ++i)
     {
     am->states[i].scores = allScores;
     allScores += rowSize;
     am->states[i].lastScores = allScores;
     allScores += rowSize;
     }
 return am;
 }
 
 void phmmMatrixFree(struct phmmMatrix **pAm)
 /* Free up memory required for an phmmMatrix and make sure
  * nobody reuses it. */
 {
 struct phmmMatrix *am = *pAm;
 if (am != NULL)
     {
     freeMem(am->states);
     freeMem(am->allCells);
     freeMem(am->allScores);
     freez(pAm);
     }
 }
 
 struct phmmState *phmmNameState(struct phmmMatrix *am, int stateIx, 
 	char *name, char emitLetter)
 /* Give a name to a state and return a pointer to it. */
 {
 struct phmmState *state;
 assert(stateIx < am->stateCount);
 state = &am->states[stateIx];
 state->name = name;
 state->emitLetter = emitLetter;
 return state;
 }
 
 static void phmmFindMatrixIx(struct phmmMatrix *am, struct phmmMommy *cell, 
 	int *retStateIx, int *retQix, int *retTix)
 /* Given a cell in matrix return state, query, and target index. */
 {
 int cellIx = cell - am->allCells;
 *retStateIx = cellIx/am->stateSize;
 cellIx %= am->stateSize;
 *retTix = cellIx / am->qDim;
 *retQix = cellIx % am->qDim;
 }
 
 static struct phmmMommy *phmmFindMommy(struct phmmMatrix *am, 
 	struct phmmMommy *baby)
 /* Find baby's mommy and return it. */
 {
 int momStateIx, qOff, tOff;
 int babyStateIx, qIx, tIx;
 UBYTE mommy;
 
 if ((mommy = baby->mommy) == phmmNullMommy)
     return NULL;
 phmmUnpackMommy(mommy, &momStateIx, &qOff, &tOff);
 phmmFindMatrixIx(am, baby, &babyStateIx, &qIx, &tIx);
 return am->states[momStateIx].cells + (tOff + tIx) * am->qDim + (qOff + qIx);
 }
 
 struct phmmAliPair
 /* Each position in alignment gets one of these. */
     {
     struct phmmAliPair *next;
     int queryIx;
     int targetIx;
     UBYTE hiddenIx;  /* If I code one of these with more than 255 states shoot me! */
     char querySym;
     char targetSym;
     char hiddenSym;
     };
 
 struct phmmAliPair *phmmTraceBack(struct phmmMatrix *am, struct phmmMommy *end)
 /* Create list of alignment pair by tracing back through matrix from end
  * state back to a start.*/
 {
 struct phmmAliPair *pairList = NULL, *pair;
 struct phmmMommy *cell, *parent = end;
 int parentSix, parentTix, parentQix;
 int sIx, tIx, qIx;
 
 phmmFindMatrixIx(am, parent, &parentSix, &parentQix, &parentTix);
 for (;;)
     {
     cell = parent;
     sIx = parentSix;
     tIx = parentTix;
     qIx = parentQix;
 
     if ((parent = phmmFindMommy(am, cell)) == NULL)
         break;
     phmmFindMatrixIx(am, parent, &parentSix, &parentQix, &parentTix);
     
     cell->mommy |= phmmMommyTraceBit;
 
     AllocVar(pair);
     pair->hiddenIx = (UBYTE)sIx;
     pair->hiddenSym = am->states[sIx].emitLetter; 
 
     if (parentQix == qIx - 1 && parentTix == tIx - 1 )
         {
         pair->queryIx = qIx-1;
         pair->querySym = am->query[qIx-1];
         pair->targetIx = tIx - 1;
         pair->targetSym = am->target[tIx-1];
         }
     else if (parentQix == qIx) /* && parentTix == tIx-1 */
         {
         pair->queryIx = -1;
         pair->querySym = '-';
         pair->targetIx = tIx-1;
         pair->targetSym = am->target[tIx-1];
         }
     else  /* parentTix == tIx  && parentQix == qIx-1 */
         {
         pair->queryIx = qIx-1;
         pair->querySym = am->query[qIx-1];
         pair->targetIx = -1;
         pair->targetSym = '-';
         }
     slAddHead(&pairList, pair);
     }
 return pairList;
 }
 
 void phmmPrintTrace(struct phmmMatrix *am, struct phmmAliPair *pairList, 
 	boolean showStates, FILE *f, boolean extraAtEnds)
 /* Print out trace to file. */
 {
 struct phmmAliPair *pair;
 #define lineLen 50
 char asyms[lineLen+1];
 char qsyms[lineLen+1];
 char tsyms[lineLen+1];
 char hsyms[lineLen+1];
 int lineIx = 0;
 int qs, ts;
 boolean gotQs = FALSE;
 
 if ((pair = pairList) == NULL)
     {
     fprintf(f, "Empty pair list\n");
     return;
     }
 
 qs = pair->queryIx;
 ts = pair->targetIx;
 
 /* Print out up to 25 bp of initial non-matching parts. */
 if (extraAtEnds)
     {
     int qStart = qs;
     int tStart = ts;
     int i;
 
     for (i= -25; i < 0; ++i)
         {
         int qIx = qStart + i;
         int tIx = tStart + i;
         qsyms[lineIx] = (qIx >= 0 ? am->query[qIx] : ' ');
         tsyms[lineIx] = (tIx >= 0 ? am->target[tIx] : ' ' );
         asyms[lineIx] = hsyms[lineIx] = ' ';
         if (qIx >= 0 || tIx >= 0)
             {
             ++lineIx;
             if (!gotQs)
                 {
                 gotQs = TRUE;
                 qs = qIx;
                 ts = tIx;
                 }
             }
         }
     }
 
 /* Print out matching parts */
 for (pair = pairList; pair != NULL; pair = pair->next)
     {
     qsyms[lineIx] = pair->querySym;
     tsyms[lineIx] = pair->targetSym;
     asyms[lineIx] = (pair->querySym == pair->targetSym ? '|' : ' ');
     hsyms[lineIx] = pair->hiddenSym;
     if (++lineIx == lineLen)
         {
         qsyms[lineIx] = asyms[lineIx] = tsyms[lineIx] = hsyms[lineIx] = 0;
         fprintf(f, "%5d %s\n", qs, qsyms);
         fprintf(f, "%5s %s\n", "", asyms);
         fprintf(f, "%5d %s\n", ts, tsyms);
         if (showStates)
             fprintf(f, "%5s %s\n", "", hsyms);
         fputc('\n', f);
         lineIx = 0;
         if (pair->next)
             {
             qs = pair->next->queryIx;
             ts = pair->next->targetIx;
             }
         }
     }
 
 /* Print out last bit - first get last pair. */
 if (extraAtEnds)
     {
     for (pair = pairList; pair->next != NULL; pair = pair->next)
         ;
 
         {
         int qIx, tIx, i;
         int residue = lineLen - lineIx;
         for (i=1; i<=residue; i++)
             {
             if ((qIx = pair->queryIx+i) >= am->querySize)
                 qsyms[lineIx] = ' ';
             else
                 qsyms[lineIx] = am->query[qIx];
             if ((tIx = pair->targetIx+i) >= am->targetSize)
                 tsyms[lineIx] = ' ';
             else
                 tsyms[lineIx] = am->target[tIx];
             asyms[lineIx] = ' ';
             hsyms[lineIx] = ' ';
             ++lineIx;
             assert(lineIx <= lineLen);
             }
         }
     }
 
 /* Print out anything not printed out yet. */
 if (lineIx != 0)
     {
     qsyms[lineIx] = asyms[lineIx] = tsyms[lineIx] = hsyms[lineIx] = 0;
     fprintf(f, "%5d %s\n", qs, qsyms);
     fprintf(f, "%5s %s\n", "", asyms);
     fprintf(f, "%5d %s\n", ts, tsyms);
     if (showStates)
         fprintf(f, "%5s %s\n", "", hsyms);
     fputc('\n', f);
     lineIx = 0;
     }
 #undef lineLen
 }
 
 struct axt *phhmTraceToAxt(struct phmmMatrix *am, struct phmmAliPair *pairList, 
 	int score, char *qName, char *tName)
 /* Convert alignment from traceback format to axt. */
 {
 struct axt *axt;
 struct phmmAliPair *pair;
 int qEnd, tEnd;
 char *qSym, *tSym;
 int i;
 
 /* Make sure got something real. */
 if ((pair = pairList) == NULL)
     return NULL;
 
 /* Allocate memory for axt. */
 AllocVar(axt);
 axt->symCount = slCount(pairList);
 axt->qSym = AllocArray(qSym, axt->symCount+1);
 axt->tSym = AllocArray(tSym, axt->symCount+1);
 
 /* Fill in basic fields. */
 axt->qName = cloneString(qName);
 axt->tName = cloneString(tName);
 axt->qStrand = axt->tStrand = '+';
 axt->qStart = qEnd = pair->queryIx;
 axt->tStart = tEnd = pair->targetIx;
 axt->score = score;
 
 /* Store alignment symbols and keep track of last symbol used. */
 for (i=0, pair = pairList; pair != NULL; pair = pair->next, ++i)
     {
     qSym[i] = pair->querySym;
     tSym[i] = pair->targetSym;
     qEnd = pair->queryIx;
     tEnd = pair->targetIx;
     }
 
 /* Store end and return. */
 axt->qEnd = qEnd + 1;
 axt->tEnd = tEnd + 1;
 return axt;
 }