gme.cpp

ncbi源码

      testEdge(e,v,A); 
      /*testEdge tests weight of meTree if loop variable 
	e is the meEdge split, places this weight in e->totalweight field */
      if (e->totalweight < w_min)
	{
	  e_min = e;
	  w_min = e->totalweight;
	}
      e = topFirstTraverse(T,e);
    }
  /*e_min now points at the meEdge we want to split*/
  GMEsplitEdge(T,v,e_min,A);
  return(T);
}

void updateSubTreeAverages(double **A, meEdge *e, meNode *v, int direction);

/*the ME version of updateAveragesMatrix does not update the entire matrix
  A, but updates A[v->index][w->index] whenever this represents an average
  of 1-distant or 2-distant subtrees*/

void GMEupdateAveragesMatrix(double **A, meEdge *e, meNode *v, meNode *newNode)
{
  meEdge *sib, *par, *left, *right;
  sib = siblingEdge(e);
  left = e->head->leftEdge;
  right = e->head->rightEdge;
  par = e->tail->parentEdge;

  /*we need to update the matrix A so all 1-distant, 2-distant, and
    3-distant averages are correct*/
  
  /*first, initialize the newNode entries*/
  /*1-distant*/
  A[newNode->index][newNode->index] = 
    (e->bottomsize*A[e->head->index][e->head->index]
     + A[v->index][e->head->index])
    / (e->bottomsize + 1);
  /*1-distant for v*/
  A[v->index][v->index] = 
    (e->bottomsize*A[e->head->index][v->index] 
     + e->topsize*A[v->index][e->head->index])
    / (e->bottomsize + e->topsize);

  /*2-distant for v,newNode*/
  A[v->index][newNode->index] = A[newNode->index][v->index] = 
    A[v->index][e->head->index];
  
  /*second 2-distant for newNode*/
  A[newNode->index][e->tail->index] = A[e->tail->index][newNode->index]
    = (e->bottomsize*A[e->head->index][e->tail->index]
       + A[v->index][e->tail->index])/(e->bottomsize + 1);
  /*third 2-distant for newNode*/
  A[newNode->index][e->head->index] = A[e->head->index][newNode->index]
    = A[e->head->index][e->head->index];
   
  if (NULL != sib)
    {
      /*fourth and last 2-distant for newNode*/
      A[newNode->index][sib->head->index] =
	A[sib->head->index][newNode->index] = 
	(e->bottomsize*A[sib->head->index][e->head->index]
	 + A[sib->head->index][v->index]) / (e->bottomsize + 1);
      updateSubTreeAverages(A,sib,v,SKEW); /*updates sib and below*/
    }
  if (NULL != par)
    {
      if (e->tail->leftEdge == e)
	updateSubTreeAverages(A,par,v,LEFT); /*updates par and above*/
      else
	updateSubTreeAverages(A,par,v,RIGHT);
    }
  if (NULL != left)
    updateSubTreeAverages(A,left,v,UP); /*updates left and below*/
  if (NULL != right)
    updateSubTreeAverages(A,right,v,UP); /*updates right and below*/  

  /*1-dist for e->head*/
  A[e->head->index][e->head->index] = 
    (e->topsize*A[e->head->index][e->head->index]
     + A[e->head->index][v->index]) / (e->topsize+1);
  /*2-dist for e->head (v,newNode,left,right)
    taken care of elsewhere*/
  /*3-dist with e->head either taken care of elsewhere (below)
    or unchanged (sib,e->tail)*/
  
  /*symmetrize the matrix (at least for distant-2 subtrees) */
  A[v->index][e->head->index] = A[e->head->index][v->index];
  /*and distant-3 subtrees*/
  A[e->tail->index][v->index] = A[v->index][e->tail->index];
  if (NULL != left)
    A[v->index][left->head->index] = A[left->head->index][v->index];
  if (NULL != right)
    A[v->index][right->head->index] = A[right->head->index][v->index];
  if (NULL != sib)
    A[v->index][sib->head->index] = A[sib->head->index][v->index];

}      
  
void GMEsplitEdge(meTree *T, meNode *v, meEdge *e, double **A)
{
  char nodelabel[NODE_LABEL_LENGTH];
  char edgelabel[EDGE_LABEL_LENGTH];
  meEdge *newPendantEdge;
  meEdge *newInternalEdge;
  meNode *newNode;
    
  sprintf(nodelabel,"I%d",T->size+1); 
  newNode = makeNewNode(nodelabel,T->size + 1);  
  
  sprintf(edgelabel,"E%d",T->size); 
  newPendantEdge = makeEdge(edgelabel,newNode,v,0.0);   
  
  sprintf(edgelabel,"E%d",T->size+1); 
  newInternalEdge = makeEdge(edgelabel,newNode,e->head,0.0);   
  
  if (verbose)
    printf("Inserting meNode %s on meEdge %s between nodes %s and %s.\n",
	   v->label,e->label,e->tail->label,e->head->label); 
  /*update the matrix of average distances*/
  /*also updates the bottomsize, topsize fields*/
  
  GMEupdateAveragesMatrix(A,e,v,newNode);

  newNode->parentEdge = e;
  e->head->parentEdge = newInternalEdge;
  v->parentEdge = newPendantEdge;
  e->head = newNode;
  
  T->size = T->size + 2;

  if (e->tail->leftEdge == e) 
    {
      newNode->leftEdge = newInternalEdge;
      newNode->rightEdge = newPendantEdge;
    }
  else
    {
      newNode->leftEdge = newInternalEdge;
      newNode->rightEdge = newPendantEdge;
    }
  
  /*assign proper topsize, bottomsize values to the two new Edges*/
  
  newPendantEdge->bottomsize = 1; 
  newPendantEdge->topsize = e->bottomsize + e->topsize;
  
  newInternalEdge->bottomsize = e->bottomsize;
  newInternalEdge->topsize = e->topsize;  /*off by one, but we adjust
					    that below*/
  
  /*and increment these fields for all other edges*/
  updateSizes(newInternalEdge,UP);
  updateSizes(e,DOWN);
}

void updateSubTreeAverages(double **A, meEdge *e, meNode *v, int direction)
     /*the monster function of this program*/
{
  meEdge *sib, *left, *right, *par;
  left = e->head->leftEdge;
  right = e->head->rightEdge;
  sib = siblingEdge(e);
  par = e->tail->parentEdge;
  switch(direction)
    {
      /*want to preserve correctness of 
	all 1-distant, 2-distant, and 3-distant averages*/	
      /*1-distant updated at meEdge splitting the two trees*/
      /*2-distant updated:
	(left->head,right->head) and (head,tail) updated at
	a given edge.  Note, NOT updating (head,sib->head)!
	(That would lead to multiple updating).*/
      /*3-distant updated: at meEdge in center of quartet*/
    case UP: /*point of insertion is above e*/
      /*1-distant average of nodes below e to 
       nodes above e*/
      A[e->head->index][e->head->index] = 
	(e->topsize*A[e->head->index][e->head->index] + 
	 A[e->head->index][v->index])/(e->topsize + 1);      
      /*2-distant average of nodes below e to 
	nodes above parent of e*/
      A[e->head->index][par->head->index] = 
	A[par->head->index][e->head->index] = 
	(par->topsize*A[par->head->index][e->head->index]
	 + A[e->head->index][v->index]) / (par->topsize + 1);
      /*must do both 3-distant averages involving par*/
      if (NULL != left)
	{
	  updateSubTreeAverages(A, left, v, UP); /*and recursive call*/
	  /*3-distant average*/
	  A[par->head->index][left->head->index]
	    = A[left->head->index][par->head->index]
	    = (par->topsize*A[par->head->index][left->head->index]
	       + A[left->head->index][v->index]) / (par->topsize + 1);
	}
      if (NULL != right)
	{
	  updateSubTreeAverages(A, right, v, UP);
	  A[par->head->index][right->head->index]
	    = A[right->head->index][par->head->index]
	    = (par->topsize*A[par->head->index][right->head->index]
	       + A[right->head->index][v->index]) / (par->topsize + 1);
	}
      break;
    case SKEW: /*point of insertion is skew to e*/
      /*1-distant average of nodes below e to 
	nodes above e*/
      A[e->head->index][e->head->index] = 
	(e->topsize*A[e->head->index][e->head->index] + 
	 A[e->head->index][v->index])/(e->topsize + 1);      
      /*no 2-distant averages to update in this case*/
      /*updating 3-distant averages involving sib*/
      if (NULL != left)
	{
	  updateSubTreeAverages(A, left, v, UP);
	  A[sib->head->index][left->head->index]
	    = A[left->head->index][sib->head->index]
	    = (sib->bottomsize*A[sib->head->index][left->head->index]
	       + A[left->head->index][v->index]) / (sib->bottomsize + 1);
	}
      if (NULL != right)
	{
	  updateSubTreeAverages(A, right, v, UP);
	  A[sib->head->index][right->head->index]
	    = A[right->head->index][sib->head->index]
	    = (sib->bottomsize*A[par->head->index][right->head->index]
	       + A[right->head->index][v->index]) / (sib->bottomsize + 1);
	}
      break;


    case LEFT: /*point of insertion is below the meEdge left*/
      /*1-distant average*/
      A[e->head->index][e->head->index] = 
	(e->bottomsize*A[e->head->index][e->head->index] + 
	 A[v->index][e->head->index])/(e->bottomsize + 1);        
      /*2-distant averages*/
      A[e->head->index][e->tail->index] = 
	A[e->tail->index][e->head->index] = 
	(e->bottomsize*A[e->head->index][e->tail->index] + 
	 A[v->index][e->tail->index])/(e->bottomsize + 1);  
      A[left->head->index][right->head->index] = 
	A[right->head->index][left->head->index] = 
	(left->bottomsize*A[right->head->index][left->head->index]
	 + A[right->head->index][v->index]) / (left->bottomsize+1);
      /*3-distant avereages involving left*/
      if (NULL != sib)
	{
	  updateSubTreeAverages(A, sib, v, SKEW);
	  A[left->head->index][sib->head->index]
	    = A[sib->head->index][left->head->index]
	    = (left->bottomsize*A[left->head->index][sib->head->index]
	       + A[sib->head->index][v->index]) / (left->bottomsize + 1);
	}
      if (NULL != par)
	{
	  if (e->tail->leftEdge == e)
	    updateSubTreeAverages(A, par, v, LEFT);
	  else
	    updateSubTreeAverages(A, par, v, RIGHT);
	  A[left->head->index][par->head->index]
	    = A[par->head->index][left->head->index]
	    = (left->bottomsize*A[left->head->index][par->head->index]
	       + A[v->index][par->head->index]) / (left->bottomsize + 1);
	}
      break;    
    case RIGHT: /*point of insertion is below the meEdge right*/
      /*1-distant average*/
      A[e->head->index][e->head->index] = 
	(e->bottomsize*A[e->head->index][e->head->index] + 
	 A[v->index][e->head->index])/(e->bottomsize + 1);        
      /*2-distant averages*/
      A[e->head->index][e->tail->index] = 
	A[e->tail->index][e->head->index] = 
	(e->bottomsize*A[e->head->index][e->tail->index] + 
	 A[v->index][e->tail->index])/(e->bottomsize + 1);  
      A[left->head->index][right->head->index] = 
	A[right->head->index][left->head->index] = 
	(right->bottomsize*A[right->head->index][left->head->index]
	 + A[left->head->index][v->index]) / (right->bottomsize+1);
      /*3-distant avereages involving right*/
      if (NULL != sib)
	{
	  updateSubTreeAverages(A, sib, v, SKEW);
	  A[right->head->index][sib->head->index]
	    = A[sib->head->index][right->head->index]
	    = (right->bottomsize*A[right->head->index][sib->head->index]
	       + A[sib->head->index][v->index]) / (right->bottomsize + 1);
	}
      if (NULL != par)
	{
	  if (e->tail->leftEdge == e)
	    updateSubTreeAverages(A, par, v, LEFT);
	  else
	    updateSubTreeAverages(A, par, v, RIGHT);
	  A[right->head->index][par->head->index]
	    = A[par->head->index][right->head->index]
	    = (right->bottomsize*A[right->head->index][par->head->index]
	       + A[v->index][par->head->index]) / (right->bottomsize + 1);
	}

      break;
    }
}

void assignBottomsize(meEdge *e)
{
  if (leaf(e->head))
    e->bottomsize = 1;
  else
    {
      assignBottomsize(e->head->leftEdge);
      assignBottomsize(e->head->rightEdge);
      e->bottomsize = e->head->leftEdge->bottomsize
	+ e->head->rightEdge->bottomsize;
    }
}

void assignTopsize(meEdge *e, int numLeaves)
{
  if (NULL != e)
    {
      e->topsize = numLeaves - e->bottomsize;
      assignTopsize(e->head->leftEdge,numLeaves);
      assignTopsize(e->head->rightEdge,numLeaves);
    }
}

void assignAllSizeFields(meTree *T)
{
  assignBottomsize(T->root->leftEdge);
  assignTopsize(T->root->leftEdge,T->size/2 + 1);
}


END_SCOPE(fastme)
END_NCBI_SCOPE

/*
 * ===========================================================================
 * $Log: gme.cpp,v $
 * Revision 1000.1  2004/06/01 18:09:56  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
 *
 * Revision 1.2  2004/05/21 21:41:03  gorelenk
 * Added PCH ncbi_pch.hpp
 *
 * Revision 1.1  2004/02/10 15:16:02  jcherry
 * Initial version
 *
 * ===========================================================================
 */